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METHOD AND SYSTEM FOR TRACKTNG MULTIPLE REGIONAL 
OBJECTS BY MULTI -DIMENSIONAL RELAXATION 



[FIELD OF THE INVENTION 
The invention relates generally to computerized 
techniques for processing data obtained from radiation 
reflections used to track multiple discrete object. 



1 CROSS-REFERENCE TO RELATED APPLICATIONS 

Thi s appl icat ion is a continuation of U.S. Patent 
Application Serial No. 09/312,036 filed May 14, 1999 which is 
a continuation of U.S. Patent Application Serial No. 
08/682,904 filed July 16, 1996, now U.S. Patent 5,959,574 
which is a continuation-in-part of U.S. Patent Application 
Serial No. 08/404,024, filed March 14, 1995 and issued July 
16, 1996 as U.S. Patent No. 5,537,119, which is a 
continuation-in-part of U.S. Patent Application Serial No. 
08/171,327 filed December 21, 1993, now U.S. Patent No. 
'"■5,406,289. 

BACKGROUND OF THE INVENTION 
a . Field of the Invention 

The invention relates generally to computerized 
techniques for processing data obtained from radar to track 
multiple discrete objects. 
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b . Description of the Background 

There are many situations where the courses of multiple 
objects in a region must be tracked. Typically, radar is used 
to scan the region and generate discrete images or "snapshots" 
5 based on sets of returns or observations. In some types of 
tracking systems, all the returns from any one object are 
represented in an image as a single point unrelated to the 
shape or size of the objects. "Tracking" is the process of 
identifying a sequence of points from a respective sequence of 

10 the images that represents the motion of an object. The 
tracking problem is difficult when there are multiple closely 
spaced objects because the objects can change [ their] speed 
and direction rapidly and move into and out of the line of 
sight for other objects. The problem is exacerbated because 

15 each set of returns may result from noise as well as echoes 
from the actual objects. The returns resulting from the noise 
are also called false positives. Likewise, the radar will not 
detect all echoes from the actual objects and this phenomena 
is called a false negative or "missed detect" error. For 

2 0 tracking airborne objects, a large distance between the radar 
and the objects diminishes the signal to noise ratio so the 
number of false positives and false negatives can be high. 
For robotic applications, the power of the radar is low and as 
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a result, the signal to noise ratio can also be low and the 
number of false positives and false negatives high. 

In view of the proximity of the objects to one another, 
varied motion of the objects and false positives and false 
5 negatives, multiple sequential images should be analyzed 
collectively to obtain enough information to properly assign 
the points to the proper tracks. Naturally, the larger the 
number of images that are analyzed, the greater the amount of 
information that must be processed. 

10 While identifying the track of an object, a kinematic 

model may be generated describing the [ object's] location, 
velocity and acceleration [may be generated! of the object . 
Such a model provides the means by which the [object's ] future 
motion of the object can be predicted. Based upon such a 

15 prediction, appropriate action may be initiated. For example, 
in a military application there is a need to track multiple 
enemy aircraft or missiles in a region to predict their 
objective, plan responses and intercept them. Alternatively, 
in a commercial air traffic control application there is a 

2 0 need to track multiple commercial aircraft around an airport 
to predict their future courses and avoid collision. 
Further, [ in] these and other applications, such as robotic 
applications, may use radar, sonar, infrared or other object 
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detecting radiation bandwidths for tracking objects. 



In 



particular, in robotic applications reflected radiation can be 
used to track a single object which moves relative to the 
robot (or vice versa) so the robot can work on the object. 

Consider the very simple example of two objects being 
tracked and no false positives or false negatives. The radar, 
after scanning at time t lt reports objects at two locations in 
a first observation set. That is, I"it1 the radar returns a set 
of two observations {o llf o 12 } . At time t 2 Fitl the radar 
returns a similar set of two observations {o 21 , o 22 ) from a 
second observation set. Suppose from prior processing that 
track data for two tracks T x and T 2 includes the locations at 
t 0 of two objects. Track T x may be extended through the points 
in the two sets of observations in any of four ways, as may 
track T 2 . The possible extensions of T 2 can be described as: 
{T lf o llf o 2l ) , {T lf o llf o 22 ) , {T lf o 12 , o 21 ) and {T lf o 12f o 22 } . 
Tracks can likewise be extended from T 2 in four possible ways^ 
including [,] {T 2 , o 12 , o 21 ) . Figure 1 illustrates these five 
(out of eight) possible tracks (to simplify the problem for 
purposes of explanation) . The five track extensions are 
labeled h llf h 12 , h 13 , h 14f and h 21 wherein h 11 is derived from {T x , 
°2i) / h 12 is derived from {T lf o llf o 22 ) , h 13 is derived from 
{T le o 12 , o 21 ) , h 14 is derived from {T lf o 12 , o 22 ) , and h 21 is 
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derived from {T 2f o llf o 21 ) . The problem of determining which 
such track extensions are the most likely or optimal is 
hereinafter known as the assignment problem. 

It is known from prior art to determine a figure of merit 
5 or cost for assigning each of the points in the images to a 
track. The figure of merit or cost is based on the likelihood 
that the point is actually part of the track. For example, 
the figure of merit or cost may be based on the distance from 
the point to an extrapolation of the track. Figure 1 

10 illustrates costs 5 21 and 5 22 for hypothetical extension h 21 and 
modeled target characteristics. The function to calculate the 
cost will normally incorporate detailed characteristics of the 
sensor, such as probability of measurement error, and track 
characteristics, such as likelihood of track maneuver. 

15 Figure 2 illustrates a two by two by two matrix, c, that 

contains the costs for each point in relation to each possible 
track. The cost matrix is indexed along one axis by the track 
number, along another axis by the image number and along the 
third axis by a point number. Thus, each position in the cost 

20 matrix lists the cost for a unique combination of points and 
a track, one point from each image. Figure 2 also illustrates 
a {0, 1} assignment matrix, z, which is defined with the same 
dimensions as the cost matrix. Setting a position in the 



5 



NUMERI .01USC2 

assignment matrix to "one" means that the equivalent position 
in the cost matrix is selected into the solution. The 
illustrated solution matrix selects the {h 14 , h 21 ) solution 
previously described. Note that for the above example of two 
5 tracks and two snapshots, the resulting cost and assignment 
matrices are three [ ] ^dimensional . As used in this patent 
application, the term "dimension" means the number of axes in 
the cost or assignment matrix while size refers to the number 
of elements along a typical axis. The costs and assignments 

10 have been grouped in matrices to facilitate computation. 

A solution to the assignment problem satisfies two 
constraints - first, the sum of the associated costs for 
assigning points to a track extension is minimized and, 
second, if no false positives or false negatives exist, then 

15 each point is assigned to one and only one track. 

When false positives exist, however, additional 
hypothetical track extensions incorporating the false 
positives will be generated. Further note that the random 
locations of false positives will, in general, not fit well 

20 with true data and such additional hypothetical track 
extensions will result in higher costs. Also note that when 
false negative errors exist, then the size of the cost matrix 
must grow to include hypothetical track extensions formulated 
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with "gaps" (i.e., data omissions where there should be 
legitimate observation data) for the false negatives. Thus, 
the second criteria must be weakened to reflect false 
positives not being as signed and also to permit the gap 
5 filler to be multiply assigned. With hypothetical cost 
calculated in this manner then the foregoing criteria for 
minimization will tend to materialize the false negatives and 
avoid the false positives. 

For a [3 -dimensional] three-dimensional problem, as is 
10 illustrated in Figure 1, but with N x (initial) tracks, N 2 
observations in scan 1, N 3 observations in scan 2, false 
positives and negatives assumed, the assignment problem can be 
formulated as: 

A A A 

(a) Minimize: [ ] ^ 1, c i^i* i^i, 
15 (b) Subject to: [ ] J2 X, z...=l, i 1 = 1 N, , (0.1) 



J2 z i si, i 2 = l, . . .N 2 . 

i 3 =i 123 



N, N, 



(d) 



[]£ E z iii sl ' ± 3 = it.-.N^ 
i^i 111213 

(e) dzni e {0 ' 1} "V z iii > 

20 [[1.0]] 
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where "c" is the cost and "z" is a point or observation 
assignment, as in Figure 2. 

The minimization equation or equivalently objective 
function I'M (0 . 1 1" . 0] ] ) (a) specifies the sum of the element by 
5 element product of the c and z matrices. The summation 
includes hypothesis representations z ±1±2± with observation 
number zero being the gap filler observation. Equation 
r fl (0 . 1 T ■ Oil ) (b) requires that each of the tracks T lf ...,T Nl be 
extended by one and only one hypothesis. Equation 

10 f n (0 . 1 [ . 01 1 ) (c) relates to each point or observation in the 
first observation set and requires that each such observation, 
except the gap filler, can only associate with one track butj. 
because of the "less than" condition^ it might not associate 
with any track. Equation f [1 (0 . 1 [ . 0] ] 1 (d) is like 

15 r fl (0 . 1 T . 01 1 ) (c) except that it is applicable to the second 
observation set. Equation r M (0 . 1 [ . 0] ] 1(e) requires that 
elements of the solution matrix z be limited to the zero and 
one values . 

The only known method to solve Problem Formulation 
20 r [1 (0 ■ 1 T . 011 ) exactly is a method called "Branch and Bound." 
This method provides a systematic ordering of the potential 
solutions so that solutions with a same partial solution are 
accessible via a branch of a tree describing all possible 
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solutions whereby the cost of unexamined solutions on a branch 
are incrementally developed as the cost for other solutions on 
the branch are determined. When the developing cost grows to 
exceed the previously known minimal cost (i.e., the bound) 
5 then enumeration of the tree branch terminates. Evaluation 
continues with a new branch. If evaluation of the cost of a 
particular branch completes, then that branch has lower cost 
than the previous bound so the new cost replaces the old 
bound. When all possible branches are evaluated or eliminated 
10 then the branch that had resulted in the last used bound is 
the solution. If we assume that N 2 = N 2 = N 3 = n and that 
branches typically evaluate to half there full length, then 
workload associated with "branch and bound" is proportional to 



15 evaluation. 

The Branch and Bound algorithm has been used in past 
research on the Traveling Salesman Problem. Messrs. Held and 
Karp showed that if the set of constraints was relaxed by a 
method of Lagrangian Multipliers (described in more detail 

2 0 below) then tight lower bounds could be developed in advance 
of enumerating any branch of the potential solution. By 
initiating the branch and bound algorithm with such a tight 
lower bound, significant performance improvements result in 



[] . [ ] This workload is unsuited to real time 
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that branches will typically evaluate to less than half their 
full length. 

Messrs. Frieze and Yadagar in dealing with a problem 
related to scheduling combinations of resources, as in job, 
5 worker and work site, showed that Problem Formulation 
r n (0 . 1 r . 01 1 ) applied. They further described a solution 
method based upon an extension of the Lagrangian Relaxation 
previously mentioned. The two critical extensions provided by 
Messrs. Frieze and Yadagar were: (1) an iterative procedure 

10 that permitted the lower bound on the solution to be improved 
(by "hill climbing" described below) and (2) the recognition 
that when the lower bound of the relaxed problem was 
maximized, then there existed a method to recover the solution 
of the non-relaxed problem in most cases using parameters 

15 determined at the maximum. The procedures attributed to 
Messrs. Frieze and Yadagar are only applicable to the [3- 
dimensionall three-dimensional problem posed by Problem 
Formulation [[] (0 . 1 [ ■ 0] ] 1 and where the cost matrix is fully 
populated. However, tracking multiple airborne objects 

20 usually requires solution of a much higher dimensional 
problem. 

Figures 1 and 2 illustrate an example where "look ahead" 
data from the second image improved the assignment accuracy 
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for the first image. Without the look ahead, and based only 
upon a simple nearest neighbor approach, the assignments in 
the first set would have been reversed. Problem Formulation 
[ [] (0 . 1 f . 01 1 ) and the prior art only permit looking ahead one 
5 image. In the prior art it was known that the accuracy of 
assignments will improve if the process looks further ahead, 
however no practical method to optimally incorporate look 
ahead data existed. Many real radar tracking problems involve 
hundreds of tracks, thousands of observations per observation 

10 set and matrices with dimensions in the range of 3 to 25 
including many images of look ahead. 

It was also known that the data assignment problem may be 
simplified (without reducing the dimension of the assignment 
problem) by eliminating from consideration for each track 

15 those points which, after considering estimated limits of 
speed and turning ability of the objects, could not physically 
be part of the track. One such technique, denoted hereinafter 
the "cone method, n defines a cone as a continuation of each 
previously determined track with the apex of the cone at the 

2 0 end of the previously defined track. The length of the cone 
is based on the estimated maximum speed of the object and the 
size of the arc of the cone is based on the estimated maximum 
turning ability of the object. Thus, the cone defines a 
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region outside of which no point could physically be part of 
the respective track. For any such points outside of the 
cones, an infinite number could be put in the cost matrix and 
a zero could be preassigned in the assignment matrix. It was 
5 known for the tracking problem that these elements will be 
very common in the cost and selection matrices (so these 
matrices are "sparse"). 

It was also known in the prior art that one or more 
tracks which are substantially separated geographically from 

10 the other tracks can be separated also in the assignment 
problem. This is done by examining the distances from each 
point to the various possible tracks. If the distances from 
one set of points are reasonably short only in relation to one 
track, then they are assigned to that track and not further 

15 considered with the remainder of the points. Similarly, if a 
larger group of points can only be assigned to a few tracks, 
then the group is considered a different assignment problem. 
Because the complexity of assignment problems increases 
dramatically with the number of possible tracks and the total 

20 number of points in each matrix, this partitioning of the 
group of points into a separate assignment problem and removal 
of these points from the matrices for the remaining points, 
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substantially reduces the complexity of the overall assignment 
problem . 

A previously known Multiple Hypothesis Testing (MHT) 
algorithm (see Chapter 10 of S.S. Blackman[,] . Multiple-Target 
5 Tracking with Radar Applications [ , Chapter 10,]^ Artech House^ 
Norwood^ MA, 1986^) related to formulation and scoring. The 
MHT procedure describes how to formulate the sparse set of all 
reasonable extension hypothesis (for Figure 1 the set 
{h 21 . . .h 24 } ) and how to calculate a cost of the hypothesis [T if 

10 o ljf o 2k ] based upon the previously calculated cost for 
hypothesis {T if o 2j } . The experience with the MHT algorithm, 
known in the prior art, is the basis for the assertion that 
look ahead through k sets of observations results in improved 
assignment of observations from the first set to the track. 

15 In theory, the MHT procedure uses the extendable costing 

procedure to defer assignment decision until the accumulated 
evidence supporting the assignment becomes overwhelming. When 
it makes the assignment decision^ it then eliminates all 
potential assignments invalidated by the decision in a process 

20 called "pruning the tree." In practice, the MHT hypothesis 
tree is limited to a fixed number of generations and the 
overwhelming evidence rule is replaced by a most likely and 
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feasible rule. This rule considers each track independently 
of others and is therefore a local decision rule. 

A general object of the present invention is to provide 
an efficient and accurate process for assigning each point 
5 object in a region from multiple images to a proper track and 
then taking an action based upon the assignments. 

A more specific object of the present invention is to 
provide a technique of the foregoing type which determines the 
solution of a Jc-dimensional assignment problem where n k n is 
10 greater than or equal to three. 

SUMMARY OF THE INVENTION 
The present invention relates to a method and apparatus 
for tracking objects. In particular, the present invention 
15 tracks movement or trajectories of objects by analyzing 
radiation reflected from the objects, the invention being 
especially useful for real-time tracking in noisy 
environments . 

In providing such a tracking capability, a region 
20 containing the objects is repeatedly scanned to generate a 
multiplicity of sequential images or data observation sets of 
the region. One or more points (or equivalently 

observations) , in each of the images or observation sets are 
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detected wherein each such observation either corresponds to 
an actual location of an object or is an erroneous data point 
due to noise. Subsequently, for each observation detected, 
figures of merit or costs are determined for assigning the 
5 observation to each of a plurality of previously determined 
tracks. * Afterwards, a first optimization problem is 
specified which includes: 

(a) a first objective function for relating the above 
mentioned costs to potential track extensions through the 

10 detected observations (or simply observations) ; and 

(b) a first collection of constraint sets wherein each 
constraint set includes constraints to be satisfied by 
the observations in a particular scan to which the 
constraint set is related. In general, there is a 

15 constraint for each observation of the scan, wherein the 

constraint indicates the number of track extensions to 
which the observation may belong. 

In particular, the first optimization problem is 
formulated, generated or defined as an M-dimensional 
2 0 assignment problem wherein there are M constraint sets in the 
first collection of constraint sets (i.e., there are M scans 
being examined) and the first objective function minimizes a 
total cost for assigning observations to various track 
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extensions wherein terms are included in the cost, such that 
the terms have the figures of merit or costs for hypothesized 
combinations of assignments of the observations to the 
tracks. Subsequently, the formulated M-dimensional assignment 
5 problem is solved by reducing the complexity of the problem by 
generating one or more optimization problems each having a 
lower dimension and then solving each lower dimension 
optimization problem. That is, the M-dimensional assignment 
problem is solved by solving a plurality of optimization 

10 problems each having a lower number of constraint sets. 

The reduction of the M-dimensional assignment problem to 
a lower dimensioned problem is accomplished by relaxing the 
constraints on the points of one or more scans thereby (■ 
permitting these points to be assigned to more than one track 

15 extension. In relaxing the constraints, terms having penalty 
factors are added into the objective function thereby 
increasing the total cost of an assignment when one or more 
points are assigned to more than one track. Thus, the 
reduction in complexity by this relaxation process is 

20 iteratively repeated until a sufficiently low dimension is 
attained such that the lower dimensional problem may be solved 
directly by known techniques. 
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In one embodiment of the invention, each Jc- dimensional 
assignment problem r2<k^M, 1 (2 < k < M) is iteratively reduced 
to a rk>11 (k-1) - dimensional problem until a [2- 
dimensionall two-dimensional problem is specified or 

5 formulated. Subsequently, the f2 -dimensional] two-dimensional 
problem formulated is solved directly and a "recovery" 
technique is used to iteratively recover an optimal or near- 
optimal solution to each /c-dimensional problem from a derived 
{k-1) [ dimensional] -dimensional problem fk=2l k= 2 ,3,4, . . .M. 

10 In performing each recovery step (of obtaining a solution 

to a ^-dimensional problem using a solution to a (k-1) - 
dimensional problem) ^ an auxiliary function [,] is utilized. 
In particular, to recover an optimal or near-optimal solution 
to a Jc-dimensional p roblem, an auxiliary ^function, 

15 W k I T , 1 (k=4 , 5 , . . . ,Ml, is specified and a region or domain is 
determined wherein this function is maximized, whereby values 
of the region determine the penalty factors of the (k-1) - 
dimensional problem such that another 12 -dimensional] two- 
dimensional problem can be formulated which determines a 

20 solution to the Jc-dimensional problem using the penalty 
factors of the {k-1) -dimensional problem. 

Each Auxiliary function W kt k=3 , . . . ,M-1, is a function 
of both lower dimensional problems , penalty factors^ and a 
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solution at the dimension k at which the penalized cost 
function is solved directly (typically a [2 -dimensional] two- 
dimensional problem) . Further, in determining, for auxiliary 
function W kf a peak region, a gradient of the auxiliary 
5 function is determined, and utilized to identify the peak 
region. Thus, gradients are used for each of the approximation 
functions W 3t W 4/ . . . , ¥ M _, which are sequentially formulated and 
used in determining penalty factors [penalty factors ] until 
rtm-ll M- 1 is used in determining the penalty factors for the 

10 [M-l dimensional 1 (M-l) -dimensional problem. Subsequently, 
once the M-dimensional problem is solved (using a [2- 
dimensionall two-dimensional problem to go from an (M-l) [ 
dimensional 1 -dimensional solution to an M-dimensional 
solution) , one or more of the following actions are taken 

15 based on the track assignments: sending a warning to aircraft 
or a ground or sea facility, controlling air traffic, 
controlling anti-aircraft or anti-missile equipment, taking 
evasive action, working on one of the objects. 

According to one feature of this first embodiment of the 

20 present invention, the following steps are also performed 
before the step of defining the auxiliary function. A 
preliminary auxiliary function is defined for each of the 
lower dimensional problems having a dimension equal or one 
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greater than the dimension at which the penalized cost 
function is solved directly. The preliminary auxiliary 
function is a function of lower order penalty values and a 
solution at the dimension at which the penalized cost function 
5 was solved directly. In determining a gradient of the 
preliminary auxiliary function, step in the direction of the 
gradient to identify a peak region of the preliminary 
auxiliary function and determine penalty factors at the peak 
region. Iteratively repeat the defining, gradient 

10 determining, stepping and peak determining steps to define 
auxiliary functions at successively higher dimensions until 
the auxiliary function [at 6-18 (k-1) 1 dimension (k-1) is 
determined. In an alternative second embodiment of the 
present invention, instead of reducing the dimentionality of 

15 the M-dimensional assignment problem by a single dimension at 
a time, a plurality of dimensions are relaxed simultaneously. 
This new strategy has the advantage that when the M- 
dimensional problem is relaxed directly to a [ 2- 
dimensional] two-dimensional assignment problem, then all 

20 computations may be performed precisely without utilizing an 
auxiliary function such as W k as in the first embodiment. More 
particularly, the second embodiment solves the first 
optimization problem (i.e., the [M- dimensional] M-dimensional 
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assignment problem) by specifying (i.e., creating, generating, 
formulating and/or defining) a second optimization problem. [ 
] The second optimization problem includes a second objective 
function and a second collection of constraint sets wherein: 
5 a) the second objective function is a combination of the 
first objective function and penalty factors or terms 
determined for incorporating ^M-m[ constraint! ) - 
constraint sets of the first optimization problem into 
the second objective function; 
10 b) the constraint sets of the second collection include only 
m of the constraint sets of the first collection of 
constraints , 

wherein r2^m^M-ll 2 < m < M -1 . Note that, once the second 
optimization problem has been specified or formulated, an 

15 optimal or near-optimal solution is determined and that 
solution is used in specifying (i.e, creating, generating, 
formulating and/or defining) a third optimization problem of 
SM-m^ dimensions (or equivalently constraint sets) . The third 
optimization problem is subsequently solved by decomposing it 

2 0 using the same procedure of this second embodiment as was used 
to decompose the first optimization problem above. Thus, a 
plurality of instantiations of the third optimization problem 
are specified, each successive instantiation having a lower 
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number of dimensions, until an instance of the third 
optimization problem is a two [ ] ^dimensional assignment 
problem which can be solved directly. Subsequently, whenever 
an instance of the third optimization problem is solved, the 
5 solution is used to recover a solution to the instance of the 
first optimization problem from which this instance of the 
third optimization was derived. Thus, an optimal or near- 
optimal solution to the original first optimization problem 
may be recovered through iteration of the above steps. . 

10 As mentioned above, the second embodiment of the present 

invention is especially advantageous when fm=2l m = 2 , since in 
this case all computations may be performed precisely and 
without utilizing auxiliary functions. 

Other features and benefits of the present invention will 

15 become apparent from the detailed description with the 
accompanying drawings contained hereinafter. 

BRIEF DESCRIPTION OF THE [FIGURES] DRAWINGS 
Figure 1 is a graph of images or data sets generated by 
20 a scan of a region and possible tracks within the images or 
data sets according to the prior art . 

Figure 2 illustrates cost and assignment matrices for the 
data sets of Figure 1 according to the prior art. 
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Figure 3 is a block diagram of the present invention. 

Figure 4 is a flow chart of a process according to the 
prior art for solving a [3 -dimensional] three-dimensional 
assignment problem. 
5 Figure 5 [I] is a flow chart of a process according to the 

present invention for solving a k-dimensional assignment 
problem where "Jc" is greater than or equal to 3. 

Figure 6 is a graph of various functions used to explain 
the present invention. 
10 Figure 7 is another block diagram of the present 

invention for solving the Jc-dimensional assignment problem 
where "k" is greater than or equal to three ( 3j_. 

Figure 8 is a flowchart describing the procedure for 
solving ag n-dimensional assignment problem according to the 
15 second embodiment of the invention. 

DETAILED DESCRIPTION OF THE [PREFERRED EMBODIMENTS] INVENTION 
Referring now to the other figures in detail wherein like 
reference numerals indicate like elements throughout the 
2 0 several views, Figure 3 illustrates a system generally- 
designated 100 for implementing the present invention. System 
100 comprises, for example, a radar station 102 (note sonar, 
microwave, infrared and other radiation bandwidths are also 
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contemplated) for scanning a region which may be, for example, 
an aerial region (in aerial surveillance applications) or a 
work region (in robotic applications) and generating signals 
indicating locations of objects within the region. The 
5 signals are input to a converter 104 which converts the 
signals to data points or observations in which each object 
(or false positive) is represented by a single point. The 
output of the converter is input to and readable by a computer 
106. As described in more detail below, the computer 106 

10 assigns the points to respective tracks, and then displays the 
tracks and extrapolations of the tracks on a monitor 110. 
Also, the computer 106 determines an appropriate action to 
take based on the tracks and track extensions. For example, 
in a commercial application at an airport, the computer can 

15 determine if two aircraft being tracked are on a collision 
course and if so, signal a transmitter 112 to warn each 
aircraft, or if a scheduled take-off will pose the risk of 
collision, delay the take-off. For a military application on 
a ship or base, the computer can determine subsequent 

20 coordinates of enemy aircraft and send the coordinates to an 
antiaircraft gun or missile 120 via a communication channel 
122. In a robotic application, the computer controls the 
robot to work on the proper object or portion of the object. 
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The invention generates Jc-dimensional matrices where k is 
the number of images or sets of observation data in the look 
ahead window plus one. Then, the invention formulates a k- 
dimensional assignment problem as in I'M ( 0 . 1 [ . 0] ] ) above . 
5 The ^-dimensional assignment problem is subsequently 

relaxed to a {k-1) -dimensional problem by incorporating one 
set of constraints into the objective function using a 
Lagrangian relaxation of this set. Given a solution of the 
(k-1) -dimensional problem, a feasible solution of the k- 

10 dimensional problem is then reconstructed. The (k-1) - 
dimensional problem is solved in a similar manner, and the 
process is repeated until it reaches the two-dimensional 
problem that can be solved exactly. The ideas behind the 
Lagrangian relaxation scheme are outlined next. 

15 Consider the integer programming problem 

Minimize v(z) = c T z 

Subject to: Az <r b, (0.2) 
Bz < 

2 0 Zt is an integer for i € I± 

[[1.0.1]] 

where the partitioning of the constraints is natural in some 
sense. Given a multiplier vector [u^Ol u £ 0 , the Lagrangian 
relaxation of [ [1 . ] JO . [1] ] 2X relative to the constraints 
25 fBz^dlBz £ d is defined to be 
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&(u) = Minimize 



{c T z+u T (Bz-d) } 



(0-3) 



5 



Subject to: 



[Az^bl Az <z b, 

Zi is an integer for i £ I± 

[[1.0.2]] 



If the constraint set fBz^dl Bz £ d is replaced by [Bz=d1 Bz 
d . the non^negativity constraint on u is removed. 
£?=c T z+u T (Bz-d) is the Lagrangian relative to the constraints 
[Bz<dl Bz < d , and hence the name Lagrangian relaxation. The 
10 following fact gives the relationship between the objective 
functions of the original and relaxed problems. 

FACT A.l. If z~ is an optimal solution to [[1.110. [11 12) , 
then 0(u) [<]£^vCz) for all [u>0l u > 0 . If an optimal solution 
z of [[l.]±0. [21131 is feasible for [ [1.1 10. [2] 121, then z is 
15 an optimal solution for [ [1 . 1 10 . [1] ] g± and &(u) [=] = v(z~) . 
This leads to the following algorithm: 

Algorithm. Construct a sequence of multipliers {u k }% =0 
converging to the solution i[ of Maximize {$>(u) : l"u^0l u > 0 ) 
and a corresponding sequence of feasible solutions {z k }^ 0 of 
20 [ [1-] 10. [1] ]21 as follows: 

1. Generate [in] an initial approximation u 0 . 

2 . Given u k , choose a search direction s k and a search 
distance a k so that &( Tu^+oisJ u , +a isj f>1 > ^(uJ . Update the 
multiplier estimate u k by [u k+1 =u k +a k s k ] = u k + C( kS k . 
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3. Given u^ +1 and a feasible solution (u^ +1 ) of 
[ [1 . ] JO . [2] ] 21, recover a feasible solution z k+1 (u k+1 ) of the 
integer programming problem f Tl . 1 (0. fl1l2) . 

4. Check the termination criteria. If the algorithm is 
5 not finished, set k=k+l and return to Step 2. Otherwise, 

terminate the algorithm. 

If J is an optimal solution of [ [1 . ] JO . [1] ] 21, then 
0(u k ) z &(u) <c v(z) £ Y(z k )± 
Since the optimal solution z~ = z "(u) of [ [1 . ] JO . [2] ] 21 is 

10 usually not a feasible solution of [ [1 . ] JO . [1] ] 21 for any 
choice of the multipliers, <P(u) is usually strictly less than 
v(z~) . The difference v(2;)-<P(u) is called the "duality gap, " 
which can be estimated by comparing the best relaxed and 
recovered solutions computed in the course of maximizing <P(u) . 

15 Thus, to eliminate the "<" in the constraint formulations 

(e.g., f n (0.1 r .01 1 ) (a) - r [] (O.l LP! 1 ) (d) ) that resulted from 
false positives, the constraint sets have been modified, as 
indicated above, to incorporate a gap filler in each 
observation set to account for false positives. Thus, a zero 

2 0 for an index [ , ] = i p /", l£D<kl (1 <> v <, k) in z k ilt _ ik in problem 
[ tl] j£. [1] ] 41 below indicates that the track extension 
represented by z* ± includes a false positive in the p th 
observation set. Note that this implies that a hypothesis be 
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formed incorporating an observation with k-1 gap fillers, 
e.g., z 0 ...oi p o ...o> i P *Q- Thus, the resulting generalization of 
Problem Formulation r M (0 . 1 \ . 01 1 ) without the "less than" 
complication within the constraints is the following k- 
5 Dimensional Assignment Problems in which rk^3l Jc £3 : 

[ Minimi ze : ] Minimize [ ] z2 • • • 22 c < i z i 4 

[ ] Subject to: []£•••£ z. .=1, [ 1 i * = 

1, • • .N, . 



.2^ ••• JL/ •■• Z-f Z i,...i~l' 

i 1= 0 i, tl =0 i k =0 



(0.4) 



i 1= 0 i jtl =0 



10 [ 

1 [ ] for ij = 1, . . .Wj 



7 and j =2, . . .Jc 



Mi M„. 



2^0 1^ = 0 

15 I,... 2^ 

[ ] []z ix---i* 6 {0 ' 1} ' Vi n± n=l / ... / ^ 

[ [1.1]] 

where c and z are similarly dimensioned matrices representing 
2 0 costs and hypothetical assignments. Note, in general, for 
tracking problems these matrices are sparse. 

After formulating the Jc-dimensional assignment problem as 
in [ [1] X0 . [1] ] 4J=/ the present invention solves the resulting 

27 
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problem so as to generate the outputs required by devices 110, 
112, 122 and 130. For each observation set Oi received from 
converter 104 at time t ± where i=l, . . . ,N lf the computer 
processes Oi in a batch together with the other observation 
5 sets Oi_ k+1 , . . . , Oi and the track Ti_ k , i.e., Tj is the set of all 
tracks that have been defined up to but not including Oj . 
(Note, bold type designations refer to the vector of elements 
for the indicated time, i.e., the set of all observations in 
the scan or tracks existing at the time, etc.) The result of 

10 this processing is the new set of tracks Ti_ k+1 and a set of 
cost weighted possible solutions indicating how the tracks 
might extend to the current time t ± . At time t 1+1 the batch 
process is repeated using the newest observation set and 
deleting the oldest. Thus, there is a moving window of 

15 observation sets which is shifted forward to always 
incorporate the most recent observation set. The effect is 
that input observation sets are reused for fk-ll (k-1) cycles 
and then on the observation set's [k-th]^ reuse each 
observation within the observation set is integrated into a 

20 track. 

Figure 7 illustrates various processes implemented upon 
receipt of each observation set. Except for the addition of 
the Jc-dimensional assignment solving process 300 and the 
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modification to scoring process 154 to build data structures 
suitable for process 300, all processes in Figure 7 are based 
upon prior art. The following processes 150 and 152 extend 
previously defined tracks hi.^j based on new observations. 
5 Gate formulation and output process 156 determines, for each 
of the previously defined tracks, a zone wherein the track may 
potentially extend based on limits of velocity, 
maneuverability and radar precision. One such technique to 
accomplish this is the cone method described previously. The 

10 definition of the zone is passed to gating process 150. When 
a new observation set Oi is received, the gating process 150 
will match each member observation with the zone for each 
member of the hypothetical set h^. After all input 
observations from Oi are processed^ the new hypothesis set h ± 

15 is generated by extending each track of the prior set of 
hypothetical tracks either with missed detect gap fillers 

or with all new observation elements satisfying the track's 
zone. This is a many [ J^tot J^many matching in that each 
hypothesis member can be extended to many new observations and 

20 each new observation can be used to extend many hypotheses. 
It, however, is not a full matching in that any hypothesis 
will neither be matched to all observations nor vice versa. 
It is this matching characteristic that leads to the sparse 
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matrices involved in the tracking process. Subsequently, 
gating 150 forwards the new hypothesis set h t to filtering 
process 152. Filtering process 152 determines a smooth curve 
for each member of h ± . Such a smooth curve is more likely than 
5 a sharp turn from each point straight to the next point. 
Further, the filtering process 152 removes small errors that 
may occur in generating observations. Note that in performing 
these tasks, the filtering process 152 preferably utilizes a 
minimization of a least squares test of the points in a track 

10 hypothesis or a Kalman [F] filtering approach. 

As noted above, the foregoing track extension process 
requires knowledge of a previous track. For the initial 
observations, the following gating process 158 and filtering 
process 160 determine the "previous track" based on initial 

15 observations. In determining the initial tracks, the points 
from the first observation set form the beginning points of 
all possible tracks . After observation data from the next 
observation set is received, sets of simple two [ l^oint 
straight line tracks are defined. Then, promotion, gate 

20 formulation, and output step 162 determines a zone in which 
future extensions are possible. Note that filtering step 160 
uses curve fitting techniques to smooth the track extensions 
depending upon the number of prior observations that have 
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accumulated in each hypothesis. Further note that promotion, 
gate formulation and output process 162 also determines when 
sufficient observations have accumulated to form a basis for 
promoting the track to processes 150 and 152 as described 
5 above . 

The output of each of the filtering processes 152 and 160 
is a set of hypothetical track extensions. Each such 
extension contains the hypothetical initial conditions (from 
the previous track) , the list of observations incorporated in 

10 the extension, and distance between each observation and the 
filtered track curve. Scoring process 154 determines the 
figure of merit or cost of an observation being an extension 
of a track. In one embodiment, the cost is based on the 
above-mentioned distance although the particular formula for 

15 determining the cost is not critical to the present invention. 
A preferred formula for determining the cost utilizes a 
negative log likelihood function in which the cost is the 
negative of the sum of: (a) the logs of the distances 
normalized by sensor standard deviation parameters, and (b) 

20 the log likelihoods for events related to [ : ] track initiation, 
track termination, track maneuver, false negatives and false 
positives. Note that track maneuvers are detected by 
comparing the previous track curve with the current extension. 
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Further note that some of the other events related to, for 
example, false negatives and false positives are detected by- 
analyzing the relative relationship of gap fillers in the 
hypothesis. Thus, after determining that one of these events 
5 occurred, a cost for it can be determined based upon suitable 
statistics tables and system input parameters. The negative 
log likelihood function is desirable because it permits 
effective integration of the useful components. Copies of the 
set of hypothetical track extensions which are scored are 

10 subsequently passed directly to one of the gate formulation 
and output steps 156 and 162. Note that the scoring process 
154 also arranges the actual scores in a sparse matrix based 
upon observation identifiers, and passes them to Jc-dimensional 
assignment problem solving process 300. 

15 The assignment solving process 300 is described below. 

Its output is simply the list of assignments which constitute 
the most likely solution of the problem described by Equation 
[ [1] [1] ]±L- Note that both gate formulation and output 
processes 156 and 162 use (at different times) the list of 

2 0 assignments to generate the updated track history T ± to 
eliminate or prune alternative previous hypotheses that are 
prohibited by the actual assignments in the list, and^ 
subsequent ly^ to output any required data. Also note that 
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when one of the gate formulation and output processes 156 and 
162 accomplish these tasks, the process will subsequently 
generate and forward the new set of gates for each remaining 
hypothesis and the processes will then be prepared to receive 
5 the next set of observations. In one embodiment, the loop 
described here will generate zones for a delayed set of 
observations rather than the subsequent set. This permits 
processes 156 and 162 to operate on even observation sets 
while the scoring step 154 and Jc-dimensional solving process 

10 operate on odd sets of observations, or vice versa. 

The assignment solving process 3 00 permits the present 
invention to operate with a window size of dimension Jc-1 for 
rk>3 1 some k > 3 . The upper limit on k depends only upon the 
computational power of the computer 106 and the response time 

15 constraints of system 100. The Jc-1 observation sets within 
the processing window plus the prior track history result in 
a Jc-dimensional Assignment Problem as described by Problem 
Formulation [[1]X0. [1]]4_L. The present invention solves this 
generalized problem including the processes required to 

2 0 consider false positives and negatives, and also the processes 
required to consider sparse matrix problem formulations. 

I . A FIRST EMBODIMENT OF THE 
DIMENSIONAL ASSIGNMENT SOLVER 3 00 
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In describing a first embodiment of the k- dimensional 
assignment solver 300, it is worthwhile to also discuss the 
process of Figure 4 which is used by the solver 300. Figure 
4 illustrates use of the Frieze and Yadagar process as shown 
5 in prior art for transforming a \3 -dimensional! three- 
dimensional assignment problem into a r2 -dimensional] two- 
dimensional assignment problem and then use a hill climbing 
algorithm to solve the \3 -dimensional! three-dimensional 
assignment problem. The solver 300 uses a Lagrangian 

10 Relaxation technique (well known in the art) to reduce the 
dimension of an original [k- dimensional] k- dimensional 
assignment problem ( [k>3] k > 3 ) down to a [3- 

dimensionall three-dimensional problem and then use the process 
of Figure 4 to solve the \3 -dimensional! three-dimensional 

15 problem. Further note that the Lagrangian Relaxation 
technique is also utilized by the process of Figure 4 and that 
in using this technique the requirement that each point is 
assigned to one and only one track is relaxed. Instead, an 
additional cost, which is equal to a respective Lagrangian 

20 Coefficient [ M ]u[ M ], is added to the cost or objective 
function [ [] (0 . 1 \ . 01 1 ) (a) whenever a point is assigned to more 
than one track. This additional cost can be picked to weight 
the significance of each constraint violation differently, so 
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this additional cost is represented as a vector of 
coefficients u which are correlated with respective 
observation points. Hill climbing will then develop a 
sequence of Lagrangian Coefficients sets designated (u 0 , ...,u j; 
5 u j+1/ . . . ,Up) . That correspond to an optimum solution of the [2- 
dimensionall two-dimensional assignment problem. The 
assignments at this optimum solution are then used to 
"recover" the assignment solution of the [3 -dimensional 1 three- 
dimensional assignment problem. 

10 In step 200 of Figure 4, initial values are selected for 

the u 0 coefficients. Because the Lagrangian Relaxation process 
is iterative, the initial values are not critical and are all 
initially selected as zero . In step 2 02 , these additional 
costs are applied to the objective function r M (0 . 1 \ . 01 1 ) (a) . 

15 With the addition of the costs ["Jut"]/ the goal is still to 
assign the points which minimize the total cost. This 
transforms Equation r M (0 . 1 [ . 0] ] 1(a) , written for fk=3 1 k= 3 
and altered to exclude mechanisms related to false positives 
and negatives, into objective function [ [2] . 1 [] ] 1 (a) . In 

20 the first iteration it is not necessary to consider the 
[ M 3u[ u ] matrix because all ["]u["] values are set to zero. To 
relax the requirement that each point be assigned to one and 
only one track, the constraint Equation \ M (0 . 1 f . 01 1 ) (d) is 
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deleted, thereby permitting points from the last image to be 
assigned to more than one track. Note that while any axis can 
be chosen for relaxation, observation constraints are 
preferably relaxed. The effect of this relaxation is to 
create a new problem which must have the same solution in the 
first two axes but which can have a fdif f eringl differ solution 
in the third axis. The result is constraints [ [2] (1 . 1 [1 1 ) (b- 
d) . 

i^O i 2 =0 123 J 3 123 

"EE va =1 ' [ 

i 2 =l i,=l 



(b) Subject to: 



(c) 
(d) 



E Va^* 1 ' * 

1 2 3 

[]Z. . . 6 {0,1} 

- L l i 2 i 3 



11 1 = 1, . . . ,N U 

11 2 = 1, . . . , J7 2st 

t 1 [IV z w, . 



15 [[2.1]] 
Step 2 04 then generates from the [3 -dimensional! three- 
dimensional problem described by Problem Formulation 
[ M (0 . 1 [ . 01 1 ) a new f 2 -dimensional] two-dimensional problem 
formulation which will have the same solution for the first 

2 0 two indices. As Problem Formulation [ [2] Xl. 1 [] ] 1 has no 
constraints on the 3 rd axis, any value within a particular 3 rd 
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axis can be used in a solution, but using anything other than 
the minimum value from any [3-rd] 3^ axis has the effect of 
increasing solution cost. Conceptually, the effect of step 
204 is to change the [ 3 -dimensional! three-dimensional arrays 
5 in Problem Formulation [ [2] Jl. 1 [] ]1 into [2 -dimensional! two- 
dimensional arrays as shown in Problem Formulation 

[ [2] 2 [] ] 1 and to generate the new [2 -dimensional] two^. 
dimensional matrix m ili2 defined as shown in Equation 

[[21H.3 [] 11. 



15 



Nl w, i'Win: 



10 (a) Minimize [:] □ £ £ | (c Ws u[ ^h ) 

ij-o i 2 =o v i, 



(b) Subject to: []£ £ z =1, [ ] i, = l,...,!^ 

i,=i i,=i ^ 

(1.2) 



*2 * ^3 



W, N, 



(C) 



n£ £ z i^ i > [ ]i 2 = i/.../^ 

(d) nz V2 e{0,ll [ ] []V z. iiz , 

[[2.2]] 

Min: arg minimize ] = Minimum { c 4 , - u. I i = l, . . . / N. }[ ] . 

i l i 2 Ji * 



[2.3] 
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] The cost or objective function for the reduced problem as 
defined by [ [2] Jj^ . 2 [] ] 1 (a) , if evaluated at all possible 
values of u is a surface over the domain of [ U 3 ]^. This 
surface is referred to as 0(u) and is non-smooth but provable 
5 convex (i.e., it has a single peak and several other critical 
characteristics which form terraces) . [Because of 1 Due to the 
convex characteristics of <£(u) , the results from solving 
Problem Formulation [ [2] iX. 2 [] ] 1 at any particular Uj can be 
used to generate a new set of coefficients u j+1 whose 

10 corresponding [Problem Formulation [2.2] problem solution is 
a] cost value is closer to the peak of <P(u) . The particular 
set of Lagrangian Coefficients that will generate the [2- 
dimensional] two-dimensional problem resulting in the maximum 
cost is designated Up. To recover the solution to the [3- 

15 dimensional! three-dimensional assignment problem requires 
solving the Equation [2 3X1.21 problem corresponding to u p . 

In step 206, the two [ ] ^dimensional problem is solved 
directly using a technique known to those skilled in the art 
such as Reverse Auction for the corresponding cost and 

20 solution values. This is the evaluation of one point on the 
surface or for the first iteration 0(u o ) . 

Thus, after this first iteration, the points have been 
assigned based on all "u" values being arbitrarily set to 
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zero. Because the "u" values have been arbitrarily assigned, 
it is unlikely that these assignments are correct and it is 
likely that further iterations are required to properly assign 
the points. Step 2 08 determines whether the points have been 
5 properly assigned after the first iteration by determining if 
for this set of assignments whether a different set of "u" 
values could result in a higher total cost. Thus, step 208 is 
implemented by determining the gradient of objective function 
[ [2] SX. 2 [] ] 1 (a) with respect to u d . If the gradient is 

10 substantially non-zero (greater than a predetermined limit) 
then the assignments are not at or near enough to the peak of 
the [$(u) ] surf ace <£(u) (decision 210), and the new set of 
Lagrangian Coefficients [are]ij| determined. 

Hill climbing Step 212 determines the u j+1 values based 

15 upon the Uj values, the direction resulting from protecting the 
previous gradient into the U 3 domain, and a step size. The 
solution value of the \2 -dimensional] two-dimensional problem 
is the set of coefficients that minimize the [2- 
dimensional] two-dimensional problem and the actual cost at the 

20 minimum. Those coefficients augmented by the coefficients 
stored in m ili2 permit the evaluation (but not the minimization) 
of the cost term in Problem Formulation [[2] SX.l [] ] X. These 
two cost terms are lower and upper bounds on the actual 
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minimized cost of the [3 -dimensional] three-dimensional 
problem, and the difference between them in combination with 
the gradient is used to compute the step size. 

With this new set of [" ]u["] values, steps 202-210 are 
5 repeated as a second iteration. Steps 212 and 202-210 are 
repeated until the gradient as a function of u determined in 
step 208 is less than the predetermined limit. This indicates 
that the u p values which locate the peak area of the <P(u) 
surface are determined and that the corresponding Problem 

10 Formulation [ [2 ] Xl. 2 [] ] 1 has been solved. Step 214 will 
attempt to use the assignments that resulted from this 
particular \2 -dimensional! two-dimensional assignment problem 
to recover the solution of the \3 -dimensionall three- 
dimensional assignment problem as described below. If the 

15 limit was chosen properly so that the values are close 

enough to the peak, this recovery will yield the proper set of 
assignments that rigidly satisfies the constraint that each 
point be assigned to one and only one track. However, if the 
[»]u[ n ] values are not close enough to the peak, then the 

20 limit value for decision 210 is reduced and the repetition of 
steps 212 and 202-210 is continued. 

Step 214 recovers the \3 -dimensionall three-dimensional 
assignment solution by using the assignment values determined 
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on the last iteration through step 208. Consider the [2- 
dimensional 1 two-dimensional z assignment matrix to have l's 
in the locations specified by the list L 1 =(a i/ jbi)f =I . If the 
[3 -dimensional 1 three-dimensional z matrix is specified by 
5 placing l's at the location indicated by the list L 2 ={a i , b if 
MaibJl-i then the result is a solution of Problem Formulation 
[ [2]ii.l [] 11. Let L 3 =(m a . b .)f =1 be the list formed by the third 
index. If each member of L 3 is unique then the L 2 solution 
satisfies the third constraint so it is a solution to Problem 

10 Formulation \ M (0 . 1 [ . 0] ] ) . When this is not the case, 
recovery determines the minimal substitutions required within 
list L 3 so that it plus L 2 will be a feasible solution, i.e., 
a solution which satisfies the constraints of a problem 
formulation, but which may not optimize the objective function 

15 of the problem formulation. This stage of the recovery 
process is formulated as a \2 -dimensional! two-dimensional 
assignment probleml":! as follows. Form a new cost matrix 
[c ifj ]f tj=1 where c itj =c aibifj for ^'=1...^ and the N ± term is the 
total number of cost elements in the selected row of the [3- 

20 dimensional] three-dimensional cost matrix. Attempt to solve 
this [2 -dimensional! two-dimensional problem for the minimum 
using two constraints sets. If a feasible solution is found 
then the result will have the same form as list L x . Replace 
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the first set of ind Texl ices by the indicated (a if pairs 
taken from list L x and the result will be a feasible solution 
of Problem Formulation [[l]l£. [l]]4l. If no feasible solution 
to the new \2 -dimensionall two-dimensional problem exists then 
5 further effort to locate the peak of &{u) is required. 

I . 1 . Generalization to a Mult i -Dimensional 
Assignment Solving Process 

Let M be a fixed integer and assume that M is the 

10 dimension of the initial assignment problem to be solved. 

Thus, initially, the result of the scoring step 154 is a M- 

dimensional Cost Matrix which is structured as a sparse matrix 

(i.e., only a small percentage of the entries in the cost and 

assignment matrices are filled or non-zero) . Individual cost 

15 elements represent the likelihood that a track [Tj^ as 

extended by the set of observations { ro f li = llO,- J i = 1 , ...,M- 

l}, is not valid. Because the matrix is sparse the list of 

cost elements is stored as a packed list, and then for each 

dimension of the matrix, a vector of a variable length list of 

20 pointers to the cost elements is generated and stored. This 

organization means that for a particular observation lo ? ] 0 ±j for 

the j th list in the i th vector will be a list of pointers to all 

hypotheses in which [o^ ] O ij= participates . This structure is 
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further explained in the following section dealing with 
problem partitioning. 

The objective of the assignment solving process is to 
select from the set of all possible combinations of track 
5 extensions a subset that satisfies two criteria. First, each 
point in the subset of combinations should be assigned to one 
and only one track and therefore, included in one and only one 
combination of the subset, and second, the total of the 
scoring sums for the combinations of the subset should be 
10 minimized. This yields the following M-dimensional equations 
where [k=M] k = M : 

(a) Minimize [:]^ v k (z k ) [ = ^£ . c*_ ijt z*_ ijt 



i 1= o 



15 



(b) [ 



• • • £ =!/ : ^ii = i,...^ 

[(c) V° V° 

for ij = 1, . . .Nj 
and j =2 , . . .k-1 

20 (d) i k = 1, . . .N k 

(e) for all i 2 . . . l k 



[3.1] U 
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i „ \ 2^f • • • 2^ /Li ' • • 2*r Z i....i y ~^-i 

for ij = l,...,Nj and j = 2, ... r k-1, 



i x -0 1^=0 
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and where c k is the cost matrix [c*l lt _ ik ] which is a function of 
the distance between the observed point z k and the smoothed 
track determined by the filtering step, and v k is the cost 
function. This set of equations is similar to the set 
5 presented in Problem Formulation [[1]X0. [1]]^L except that it 
includes the subscript and superscript k notation. Thus, in 
solving the M-dimensional Assignment Problem the invention 
reduces this problem to an [M-l dimensional! (M-l) -dimensional 
Assignment Problem and then to an [N-2 dimensional! (N-2) - 

10 dimensional Assignment Problem, etc. Further, the symbol 
Jce{3,...,M} customizes Problem Formulation [ [3] JjL. [1] ] £L to a 
particular relaxation level. That is, the notation is used to 
reference data from levels relatively removed as in c k+1 are the 
cost coefficients which existed prior to this level of relaxed 

15 coefficients c k . Note that actual observations are numbered 
from 1 to N ± ± where N ± is the number of observations in 
observation set i. Further note that the added [0]O 
observation in each set of observations is the unconstrained 
"gap filler." This element serves as a filler in substituting 

20 for missed detects, and a sequence of these elements including 
only one true observation represents the possibility that the 
observation is a false positive. Also note that by being 
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unconstrained a gap filler may be used in as many hypotheses 
as required. 

While direct solution to [ [3] XI. [1] ] 11 would give the 
precise assignment, the solution of /c-dimensional equations 
5 directly for large k is too complex and time consuming for 
practice. Thus, the present invention solves this problem 
indirectly. 

The following is a short description of many aspects of 
the present invention and includes some steps according to the 

10 prior art. The first step in solving the problem indirectly 
is to reduce the complexity of the problem by the previously 
known and discussed Lagrangian Relaxation technique. 
According to the Lagrangian Relaxation technique, the absolute 
requirement that each point is assigned to one and only one 

15 track is relaxed such that for some one image, points can be 
assigned to more than one track. However, a penalty based on 
a respective Lagrangian Coefficient u k is added to the cost 
function when a point in the image is assigned to more than 
one track. The Lagrangian Relaxation technique reduces the 

20 complexity or "dimension" of the formulation of the assignment 
problem because constraints on one observation set are 
relaxed. Thus, the Lagrangian Relaxation is performed 
iteratively to repeatedly reduce the dimension until a [2- 
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dimensional 1 two-dimensional penalized cost function problem 
results as in Problem Formulation [ [2] J^. 1 [] ] 1. This [2- 
dimensional] two-dimensional problem is solved then directly by 
a previously known technique such as Reverse Auction. The 
5 penalized cost function for the [2 -dimensional! two-dimensional 
problem defines a valley or convex shaped surface which is a 
function of various sets of { [u k |k=3] u k | k=3 , . . . ,M) penalty 
values and one set of assignments for the points in two 
dimensions. That is, for each particular u 3 there is a 

10 corresponding [2 -dimensional] two-dimensional penalized cost 
function problem and its solution. Note that the solution of 
the \2 -dimensional! two-dimensional penalized cost function 
problem identifies the set of assignments for the particular 
u 3 values that minimize the penalized cost function. However, 

15 these assignments are not likely to be optimum for any higher 
dimensional problem because they were based on an initial 
arbitrary set of u k values. Therefore, the next step is to 
determine the optimum assignments for the related [3- 
dimensional] three-dimensional penalized cost function problem. 

20 There exists a [2 -dimensional 1 two-dimensional hill shaped 
function 0 which is a graph of the minimums of all penalized 
cost functions at various sets of assignments in [2- 
dimensions] two dimensions . For the \3 -dimensional] three- 
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dimensional problem, the [<t> ] function <P can be defined based 
on the solution to the foregoing [2 -dimensional! two- 
dimensional penalized cost function. By using the current u 3 
values and the { fu k |k>3l u k | k > 3 ) values originally assigned[. 
5 Then], the gradient of the hill-shaped [<J> 1 function $> is 
determined, which [gradient] points toward the peak of the 
hill. By using the gradient and u 3 values previously selected 
for the one point on the hill (corresponding to the minimum of 
the penalized cost [<E» ] function <P ) as a starting point, the 

10 u 3 values can be found for which the corresponding problem will 
result in the peak of the [<5 ] function 0. The solution of the 
corresponding \2 -dimensional] two-dimensional problem is the 
proper values for [2] two of the f 3 1 three sets of indices in 
the r3 -dimensional] three-dimensional problem. These solution 

15 indices can select a subsection of the cost array which maps 
to a [2 -dimensional] two-dimensional array. The set of indices 
which minimize the \2 -dimensional] two-dimensional assignment 
problem based on that array corresponds to the proper 
assignment of points in the third dimension. The foregoing 

2 0 "recovery" process was known in the prior art, but it is 
modified here to adjust for the sparse matrix characteristic. 
The next task is to recover the solution of the proper u 4 
values for the r4-dimensional1 four-dimensional problem. The 
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foregoing hill climbing process will not work again because 
the foregoing hill climbing process when required to locate 
the f4 -dimensional l f our-dimensional u 4 values for the peak of 
<P* requires the exact definition of the [<E> 3 ] function & (as 
5 was available in the case of <t?) or an always less than 
approximation of , whereas the iteration can result in a 
greater than approximation of # 3 . According to the present 
invention, another f3 -dimensional! three-dimensional function 
is defined which is a "less than approximation" the [3- 

10 dimensional 4> 1 three-dimensional function <P and which can be 
defined based on the solution to the \2 -dimensional! two- 
dimensional penalized cost function and the previously 
assigned and determined u k values. Next, the gradient of the 
[W 1 function W is found and hill climbing technique used to 

15 determine the u 4 values at the peak. Each selected u 4 set 
results in a new \3 -dimensionall three-dimensional problem and 
requires the \2 -dimensionall two-dimensional hill climbing 
based upon new \2 -dimensionall two-dimensional problems. At 
the peak of the [3 -dimensional W 1 three-dimensional function 

20 W, the solution values are a subset of the values required for 
the f4 -dimensionall four-dimensional solution. Recovery 
processing extends the subset to a complete solution.. This 
process is repeated iteratively until the u k values that result 
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in a corresponding solution at the peak of the highest order 
[W and <& 1 functions W and <P are found. The final recovery 
process then results in the solution of Jc-dimensional problem. 
Figure 5 illustrates process 300 in more detail. 



1.2. Problem Formulation 



In problem formulation step 310, all data structures for 
the subsequent steps are allocated and linked by pointers as 
required for execution efficiency. The incoming problem is 

10 partitioned as described in the subsequent section. This 
partitioning has the effect of dividing the incoming problem 
into a set of independent problems and thus reducing the total 
workload. The partitioning process depends only on the actual 
cost matrix so the partitioning can and is performed for all 

15 levels of the relaxation process. 



1 . 2 . 1 . [ ] _Relaxation and Recovery 
Step 320 begins the Lagrangian Relaxation process for 
reducing the M-dimensional problem by selecting all Lagrangian 
20 Coefficient u M penalty values initially equal to zero. The 
Lagrangian Coefficients associated with the 14 th constraint set 
are a[ N M +1 element] n (AL+1) -element vector. The reduction of 
this M-dimensional problem in step 322 to a [M-l 
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dimensionall (M-l) -dimensional problem uses the two step 
process described above. First, a penalty based on the value 
of the respective u M coefficient is added to the cost function 
when a point is assigned to more than one track and then the 
5 resultant cost function is minimized. However, during the 
first iteration, the penalty is zero because all u M values are 
initially set to zero. Second, the requirement that no point 
from any image can be . assigned to more than one track is 
relaxed for one of the images. In the extreme this would 

10 allow a point from the relaxed image to be associate with 
every track. However, the effect of the previous penalty 
would probably mean that such an association would not 
minimize the cost. The effect of the two steps in combination 
is to remove a hard constraint while adding the penalty to the 

15 cost function so that it operates like a soft constraint. For 
step 322 this two [ ]^step process results in the following 
penalized cost function problem with k=M: 
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(a) <P k (u k ) = Minimize : (p k (u k ,z k ) 



A/ • • • /L/ Z ij. . .ijj 1 



£ • • • £ £ u i 



(b) Subject to: z2 ■■■z2 z i 1 ...i k = ^-' t 1 = 1, 



i 2 =0 



(1-5) 



10 



(c) 



(d) 



[ 



Z-f • " • Z-r Z-f • • * Z-f Z i 1 ...i Jt - 



i x -0 



[ ] for ij = 1, . . .Nj 
[ ] and j-2 , . . 

*£...i e fo,l} [ ]_for all i 2 .,.l^ 



[3.2] 

] Because the constraint on single assignment of elements from 
the last image has been eliminated, a[ M-l dimensional] n (M- 
1) -dimensional problem can be developed by eliminating some of 

15 the possible assignments. As shown in Equations 

[ [3] s il. [3 ] ] 6Jl, this is done by selecting the smallest cost 
element from each of the M th axis vectors of the cost matrix. 
Reduction in this manner yields a new, lower [ J^order 
penalized cost function defined by Equations [ [3] (1 . [3] ] 6) 

2 0 which has the same minimum cost as does the objective function 
defined by [ [3] _[1 . [2] ] 51 (a) above. 
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[ ] The cost vectors are selected as follows. [ ] Define 

[an index array m^..^.! and] a new cost array by: 

c^...^ Minimum {c*_^-u 2 * | ^ = 0,1, . . (1.6) 
for (i lf . . . , i^) ^ (0, . . . , 0) , 



c 0 



?.o = £ q min j 0 ' Co* . .o " "i* }• 



[ [3.3] 
1 The resulting [M-l dimensional] (M-l) -dimensional problem 
5 is (where k=M) : 

0 k (u k ) = Minimize: v k _ 1 (z k ' 1 )= c £ ~\ .i k x z i~^ .i k 1 

Subject to: £ . . . 2 z i" 1 i = 1 ' i 2 = 1, . . -^i, (1.7) 

ijt =o 



2L, - * * • • • z i 1 . . .I* -1 ' 



i 1= 0 i J+1 -0 i k =0 



for ij = 2, . . .Nj 

10 and j-2, . . .k-1, 



^ ... ^ z ij> - 1, . . .N k . lt 



V 0 i*- 2 -° 



z V..i t 6 * 0 ' ^ for all ^...lw 

[ [3.4] 
]Assignment Problem [ [3] Xl. [1] ] 4_)_ and Problem Formulation 
15 [ [3] J^. [4] ] 71 differ only in the dimension M vs. M-l, 
respectively. An optimal solution to [[3]_[1.- [4]]7Jl is also an 
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optimal solution to equation [[3]li. [2]] JLL- This relationship 
is the basis for an iterative sequence of reductions indicated 
by steps 320-332 through 330-332 and 200-204 in which the 
penalized cost function problem is reduced to a two [ ]^ 
5 dimensional problem. As these formula will be used to 
describe the processing at all levels, the lowercase k is used 
except where specific reference to the top level is needed. 
In step 206, the \2 -dimensional! two-dimensional penalized cost 
function is solved directly by the prior art Reverse Auction 

10 technique. Each execution of 206 produces two results, the 
set of z 2 values that minimize the problem and the cost that 
results v 2 when these z values are substituted into the 
objective function r T31 (1 . [4] 1 7) (a) ■ 

In step 208, according to the prior art, solution z 2 

15 values are substituted into the f 2 -dimensional! two-dimensional 
derivative of the [<I> 2 1 surface &<>. The result indicates how 
the value of u 3 should be adjusted so as to perform the hill 
climbing function. As was previously described the objective 
is to produce a sequence of u- values which ends when the 

20 value is in the domain of the peak of the [$ 1 surface <P. The 
section "Determining Effective Gradient" describes how new 
values are computed and how it is determined that the current 
Ui points to the peak of <P 2 . When no further adjustment is 
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required the flow moves to step 214 which will attempt to 
recover the \3 -dimensionall three-dimensional solution as 
previously described. When further adjustment is required 
then the flow progresses to step 212 and the new values of u k 
5 are computed. At the ^-dimensional! two-dimensional level the 
method of the prior art could be used for the hill climbing 
procedure. However, it is not practical to use this prior art 
hill climbing technique to determine the updated Lagrangian 
Coefficients u k or the Max on the next (or any) higher order [0 

10 ] surface <P because the next (or any) higher dimensional [$ 
1 function <P cannot be defined based on known information. 

Instead, the present invention defines a new function 
based on known information which is useful for hill climbing 
from the third to the fourth and above dimensions, i.e., u k 

15 values which result in z values that are closer to the proper 
z values for the highest Jc-dimension . This hill climbing 
process (which is different than that used in the prior art of 
Figure 4 for recovering only the three [ ] ^dimensional 
solution) is used iteratively at all lower dimensions 7c of the 

20 M-dimensional problem (including the \3 -dimensionall three- 
dimensional level where it replaces prior art) even when k is 
much larger than three. Figure 6 helps to explain this new 
hill climbing technique and illustrates the original k- 
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dimensional cost function v k of Problem Formulation 
[[3]li. [l]]^!. However, the actual £>dimensional cost surface 
v k defined by [[3]^. [1] ] 41(a) comprises scalar values at each 
point described by Jc-dimensional vectors and as such can not 
5 be drawn. Nevertheless, for illustration purposes only, 
Figure 6 ignores the reality of ordering vectors and 
illustrates a concave function v k (z k ) to represent Equation 
[3] Ji. [1] £1. The [ v k ] surface V l. is illustrated as being 
smooth to simplify the explanation although actually it can be 
10 imagined to be terraced. The goal of the assignment problem 
is to find the values of ~z k ; these values minimize the k- 
dimensional cost function v^. 

For purposes of explanation, assume that in Figure 6, A:=4 
(the procedure is used for all rk> 3 1 Jc >3 ) . This problem is 
15 reduced by two iterations of Lagrangian Relaxation to a [2- 
dimensionall two-dimensional penalized cost function . This 
cost function, and all other cost functions described below^ 
are also non- smooth and continuous but are illustrated in 
Figure 6 as smooth for explanation purposes. Solving the 
20 problem results in one set of z 2 assignments and the value of 
<P( k . 1} . at the point u-j . A series of functions 

[$_, x , ...,<|>^ ]<Pj r • • • / 0i each generated from a 

different u 3 are shown. The particular series illustrates the 
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process of locating the peak of 0( k . 1 ) i » The [2 -dimensional! two- 
dimensional penalized cost functions 

[c|>_. 1 f ... / <t)^ x ] 0j , . . . , $i can be solved directly. 
Each such solution provides the information required to 
5 calculate the next u 3 value. Each iteration of the hill 
climbing improves the selection of u 3 values, i.e., yields <p 2 
problem whose solution is closer to those at the solution of 
the 0 3 problem. The result of solving is values that are on 
both ^ > (k . 1)i and & k . Figure 6 illustrates the [$ k 1 surface 

10 which comprises the minimums of all Jc-dimensional penalized 
cost function [$ k 1 surfaces i.e., if the Problem 

Formulations [ [3] 11. [2] ] JX and [ [3] JX . [4] ] 21 were solved at 
all possible values of u k the function <t> k (u k ) would result. The 
[$ k ] surface <£i is always less than the [v k ] surface v k _ except 

15 at the peak as described in Equation [3]jLl. [5]8X and its 
maximum occurs where the [v k 1 surface w is minimum. Because 
the [<I> k 1 surface represents the minimum [s] of the [cj) k 

] surface <P U . any point on the [$> k ] surface can mapped to 
the [v k ] surface . The [ O k ] function & k provides a lower 

20 bound on the minimization problem described by 
[ [3] XL. [1] ] 4X (a) . Let ~z k be the unknown solution to Problem 
Formulation [ [3] XL • [U ] 4X and note that: 
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[[3.5]] 



[]* k (u k )&v k ([z k ]z k )z v k _ x {z k ) . 



(1.8) 



Consequently, the [z k 1 values z k at the peak of & k (i.e., the 
greatest lower bound on the cost of the relaxed problem) , can 
5 be substituted into the k- dimensional penalized cost function 



present invention attempts to find the maximum on the [<t> k 
] surface However, it is not possible to use the prior art 

hill climbing to hill climb to the peak of & k because the 

10 definition of <P k requires exact knowledge of lower order [$ 
] f unctions = g. As the solution of <ff{ is not the exact [ the] 
solution, in that higher order [u lvalues of u are not yet 
optimized, its value can be larger than the true peak of & k . 
As such it is not a lower bound on and it can not be used to 

15 recover the required solution. 

Instead, the present invention defines all auxiliary 
functions ' which fisl are based only on the solution to the 
\2 -dimensional] two-dimensional penalized cost function 
problem, lower order values z k and values [and ] u k [ values] 

20 determined previously by the reduction process. The [W k 
] function P L is a less than approximation of <P k/ and its 
gradient is used for hill climbing to its peak. The [zl^J* 
lvalues z k at the peak of the [W k 1 function ff l are then 



to determine the proper assignments. 



Consequently, the 
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substituted into Problem Formulation [ [3]li. [2] ] 5JL to 
determine the proper assignments. To define the [¥ k ] function 
W k , the present invention explicitly makes the function ^(u*) 
a function of all higher order sets of Lagrangian Coefficients 
5 with the expanded notation: <P k (u k ; u** 1 , . . . ,u k ) . Then, a new set 
of [W 1 functions W. is defined recursively, using the ^'s 
domain [ . ] z_ 

[ [3.6] (a)], 

Vu 3 ) = v 2 + z2 ui =«>3(u 3 ;uS • • -,U K ) , 



10 



(1.9) 



i,=0 



where v 2 is the solution value for the most recent [2- 
dimensionall two-dimensional penalized cost function problem. 
For [k>3]/c_>3 



W. (u 3 , . . .,u* _1 ;u*) = < 



<S> k (u k ;u k+1 , ...,u K ), 

if known, 

W.. 1 (u 3 / ...,u ,c - 2 ;u*- 1 )+2] "iV 
otherwise. 



(1.10) 



[ [3.6] (b) 

15 ] From the definitions of & k and v k (Problem Formulation 

[ [3 ] XI. [4] ] compared with Problem Formulation 
[[3]±1. [2]]5J_) : 
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<!> k (u k ;u k+1 ,...,u M ) = v k _ 1 (z k - 1 )+22 ui, 



it follows that: 

^ 3 (" 3 ) = V 2+T, Ui 3 =^ 3 (u 3 ;u 4 f . . .,u M )i v 3 3 (z 3 ) 



and with that Equation [3]Xl. [53 8^ is extended to: 
]W k (u 3 , . . .,u*-^u*)s<^(u*;u* +1 , . . u M )l v k (u k )z v k (u k ) (1.11 

[ [3.7] 
5 ]This relationship means that either <P k or W k may be used in 
hill climbing to update the Lagrangian Coefficients u. <P k is 
the preferred choice, however it is only available when the 
solution to Problem Formulation [ [3] SX> [2] ] jy. is a feasible 
solution to Problem Formulation [ [3] Xi. [1] ] 4^ (as in hill 
10 climbing from the second to third dimension which is why prior 
art worked) . For simplicity in implementation, the function 
W k is defined so that it equals <P k when either function could 
be used. It is therefore always used, even for hill climbing 
from the second dimension. The use of the 1 function W 
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which is based on previously assigned or determined u values 
and not higher order u values which are not yet determined, is 
an important feature of the present invention. 

5 1.2.2. Determining Effective Gradient 

After the V¥ k ] function is defined, the next steps of 
hill climbing/peak detection are to determine the gradient of 
the [W k ] function W kt determine an increasing portion of the 
gradient and then move up the [^F k ] surface W k in the direction 

10 of this increasing portion of the gradient. As shown in 
Figure 6, any upward step on the [W k 1 surface W^, for example^ 
to the minimum of the [0 k j]^j will yield a new set of values 
u k [ values] (to the "left") that is closer to the ideal set of 
values u k [ values] which correspond to the minimum of the [v k 

15 ] function v k . While it is possible to iteratively step up this 
[Y k 1 surface W k with steps of fixed size and then determine if 
the peak has been reached, the present invention optimizes 
this process by determining the single step size from the 
starting point at the minimum of that will jump to the 

20 peak and then [calculating] calculate the values u k [ values] at 
the peak. Once the [f."u" si values u at the peak [at ] of ffl 
are determined, then the values u k [ values] can be substituted 
into Problem Formulation [[3]Xl. [2]]|lL to determine the proper 
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assignment. (However, in a more realistic example, where k is 
much greater than three, then the values u k [ values] at the 
peak of the ] function ¥ along with the lower [ l^order 

values u k [ values] and those assigned and yielded by the 
5 reduction steps are used to define the next higher level [W 
] function ¥. This definition of a higher order ] function 
¥ and hill climbing process are repeated iteratively until 
¥ k , [ and] the peak of ¥ k ^ = and the values u k [ values] at the peak 
of ¥ k are identified.) The following is a description of how 

10 to determine the gradient of each [W ] surface ¥ and how to 
determine the single step size to jump to the peak from the 
starting point on each [W 1 surface ¥. 

As noted above, each [W 1 surface ¥ is non-smooth. 
Therefore, if a gradient [were] was taken at a single point, 

15 the gradient may not point toward the peak. Therefore, 
several gradients (a "bundle") are determined at several 
points on the [W k ] surface in the region of the starting 
point (i.e., minimum of <p k .) and then averaged. Statistically, 
the averaged result should point toward the peak of the 

20 ] surface ¥ k . Wolfe ! s Conjugate Subgradient Algorithm 

( [Wolfe75, Wolfe79] P. Wolfe. A method of conjugate 
subgradients for minimizing non-dif f erentiable functions. 
Mathematical Programming Study, 3:145-173, 1975; P. Wolfe. 
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Finding the nearest point in a polvtope. Mathematical 
Programming- Study, 11:128-149, 1976 ) for minimization was 
previously known in another environment to determine a 
gradient of a non- smooth surface using multiple subgradients 
5 and can be used with modification in the present invention. 
[Wolfed algorithm is further described in "A Method of 
Conjugate Subgradients for Minimizing Nondif f erentiable 
Functions" page 147-173 published by Mathematical Programming 
Study 3 in 1975 and "Finding the Nearest Point in a Polytope" 

10 page 12 8-149 published by Mathematical Programming Study 11 in 
1976. ] The modification to Wolfe's algorithm uses the 
information generated for ^(u k3 , . . . ,u kk ~ 2 ; u kk_1 ) as the basis for 
calculating the subgradients. The definition of a subgradient 
v of W k (u k3 , . . . , u kk ) is any member of the subdif f erential set 

15 defined as: 

[5^ k ] ^Sc(u)={ fv6R Nk+i n v E R N * +1 I (SUu k3 . . . . / u kk - 1 ;u') [-^ k ]^ 

V k (u k2 t . . . / u k " 1 ;u k ) ) [£v T ] 

_^ ^(u'-u k ) [ V u'ER Nk+1 l V u'g R N ^M . 

2 0 [(] where [v t ]yI is the transpose of vt)]^ 
[ 

1 Next . a subgradient vector is determined from this 

function. If z k is the solution of Problem Formulation 
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[ [3]ll. [2] then differentiating W k (u 3 , . . . , u* 1 ; u k ) with 

respect to u k k and evaluating the result with respect to the 
current selection matrix z k yields a subgradient vector: 

-EH- E ...S[UM] urn [ Zi i Ji k =l, ...,N k ))]z* ' , \ 



5 The iterative nature of the solution process at each dimension 
yields a set of such subgradients . Except for a situation 
described below^ where the resultant averaged gradient does 
not point toward the peak, the most recent set of such 
subgradients are saved and used as the "bundle" for the peak 

10 finding process for this dimension. For example, at the [k] 
level k there is a bundle of subgradients of the [¥ k ] surface 
JPJ^ near the minimum of the [0 ki ] surface 0* de t e rmi ned as a 
result of solving Problem Formulation [ [3 ] J^L . [1] ] 4^ at all 
lower levels. This bundle can be averaged to approximate the 

15 gradient. Alternately, the previous bundle can be discarded 
so as to use the new value to initiate a new bundle. This 
choice provides a way to adjust the process to differing 
classes of problems, i.e., when data is being derived from two 
sensors and the relaxation proceeds from data derived from one 

2 0 sensor to the other then the prior relaxation data for the 
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first sensor could be detrimental to performance on the second 
sensors data. 

1.2.3. [ 1 _Step Size and Termination Criterion 
5 After the average gradient of the [V k ] surface ffl is 

determined, the next step is to determine a single step that 
will jump to the peak of the [*F k 1 surface ¥ k . The basic 
strategy is [to 1 first to specify an arbitrary step size, and 
then calculate the value of W k at this step size in the 

10 direction of the gradient. If the [W k 1 value of W k is larger 
than the previous one, this probably means that the step has 
not yet reached the peak. Consequently, the step size is 
doubled and a new [W k lvalue of is determined and compared 
to the previous one. This process is repeated until the new 

15 [W k lvalue of is less than the previous one. At that time, 
the step has gone too far and crossed over the peak. 
Consequently, the last doubling is rolled-back and that step 
size is used for the iteration. If the initial estimated W k 
is less than previous value then the step size is decreased by 

20 50%. If this still results in a smaller [W k lvalue of then 
the last step is rolled back and the previous step size is 
decreased by 25%. The following is a more detailed 
description of this process. 
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With a suitable bundle of subgradients determined as just 
described [.] j_ Wolfe's algorithm can be used to determine the 
effective subgradient [(]d[)] and the [Ulupgraded value u* +1 . 
From the previous iteration, or from an initial condition, 
5 there exists a step length value [(]t[)]. The value, 
[ 

] u + = u* + td 

[ 

] is calculated as an estimate of u* +1 . To determine if the 
10 current step size is valid [the] we evaluate ^(u 3 , . . . ,u k " 2 ; uj . 
If the result represents an improvement then we double the 
step size. Otherwise we halve the step size. In either case 
a new u + is calculated. The doubling or halving continues 
until the step becomes too large to improve the result, or 
15 until it becomes small enough to not degrade the result. The 
resulting suitable step size is saved with d as part of the 
subgradient bundle. The last acceptable u + is assigned to u* +1 . 

Three distinct criteri[on]a are used to determine when u* 
is close enough to u k : 
20 1. The Wolfe's algorithm criterion of d=0 given that the 
test has been repeated with the bundle containing only 
the most recent subgradient. 
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2. The difference between the lower bound ^(u k ) and the 
upper bound v^(z k , u k ) being less than a preset relative 
threshold. (Six percent was found to be an effective 
threshold for radar tracking problems.) 

3. An iteration count being exceeded. 

The use of limits on the iteration are particularly 
effective for iterations at the level 3 through Jn-1[ 
levels,]^ as these iterations will be repeated so as to 
resolve higher order coefficient sets. With limited 

iterations the process is in general robust enough to improve 
the estimate of upper order Lagrangian Coefficients. By- 
limiting the iteration counts then the total processing time 
for the algorithm becomes deterministic. That characteristic 
means the process can be effective for real time problems such 
as radar, where the results of the last scan of the 
environment must be processed prior to the next scan being 
received. 

1.2.4. [ 1 _Solution Recovery 
The following process determines if the u k values at what 
is believed to be the peak of the ] function adequately 

approximate the u k values at the peak of the corresponding [<3> k 
] function <P k . This is done by determining if the corresponding 
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penalized cost function is minimized. Thus, the u k values at 
the peak of the [*P k ] function W k are first substituted into 
Problem Formulation [ [3] Xl . [2] ] 51 to determine a set of z 
assignments for k-1 points. During the foregoing reduction 
5 process, Problem Formulation [[3]J^. [4]] 71 yielded a tentative 
list of k selected z points that can be described by their 

indices as: [{(i^ • * • i k-l)}j= - " ' ij ^ _1 )} W -i' where N o is the 

number of cost elements selected into the solution. One 
possible solution of Problem Formulation [[3]XI. [l]]4l is the 

10 solution of Problem Formulation [ [3] S±. [2] ] ELL which is 
described as [{(ij . . . ij^m ^ . ± J} ^ {( i{ . . . i^ia*. . .J}^ 
with = as it was defined in Problem Formulation [[3]Xl. \3] 1 6) . 
If this solution satisfies the [k-th]k^ constraint set^ then 
it is the optimal solution. 

15 However, if the solution for Problem Formulation 

[ [3] 11. [2] ] 51 is not feasible (decision 355), then the 
following adjustment process determines if a solution exists 
which satisfies the [k-th]^^ constraint while retaining the 
assignments made in solving Problem Formulation [[3]^!. K1 1 7) . 

20 To do this, a [2 -dimensional] two-dimensional cost matrix is 
defined based upon all observations from the [k-th]^ set 
which could be used to extend the relaxed solution. 
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1_ K 



for 1 = 0, . . . ,N k 
"Vi' 1 and j = 0, . . . , N Q . 



(1.12) 



[[3.8]] 

If the resulting [2 -dimensional] two-dimensional assignment 
problem, 



Minimi ze [ : ] [ 1 X) £ h i i 



j=0 1=0 " J 



Subject to: [l^L^^l t 1 i=l, . . . ,N n , 

1 = 0 3 

(1.13) 



1 = 0 

[1^6 {0,1} [ ] j = 0,...,N 0 , 1 = 0,..., N k 

[ [3.9] 
] has a feasible solution^ then the indices of that solution 
map to the solution of Problem Formulation [ [3]^. [1] ]4L for 
the Jc-dimensional problem. The first index in each resultant 
solution entry is the pointer back to an element of the [ 
]list. That element supplies the first k-1 indices of the 
solution. The second index of the solution to the recovery- 
problem is the k th index of the solution. Together these 
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indexes specify the values of z k that solve Problem Formulation 
[ [3] 11. [1] ]4l at the k th level. 

If Problem Formulation \ [3 . 91 1 (1 . 13) does not have a 
feasible solution^ then the value of which was thought to 
5 represent u£ is not representative of the actual peak and 
further iteration at the k** level is required. This decision 
represent the branch path from step 214 and equivalent steps. 

10 1.3. Partitioning 

A partitioning process is used to divide the cost matrix 
that results from the Scoring step 154 into as many 
independent problems as possible prior to beginning the 
relaxation solution. The partitioning process is included 

15 with the Problem Formulation Step 310. The result of 
partitioning is a set of problems to be solved, i.e., there 
will be Pj problems that consist of a single hypothesis, p 2 
problems that consist of two hypothesis, etc. Each such 
problem is a group in that one or more observations or tracks 

2 0 are shared between members of the group. 

In partitioning to groups no consideration is given to 
the actual cost values. The analysis depends strictly on the 
basis of two or more cost elements sharing the same specific 
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axis of the cost matrix. In a [2 -dimensional! two-dimensional 
case two cost elements must be in the same group if they share 
a row or a column. If the two elements are in the same row, 
then each other cost element that is also in the same row, as 
5 well as any cost elements that are in columns occupied by 
members of the row must be included in the group. The process 
continues recursively. In literature it is referred to as 
"Constructing a Spanning Forest." The 7c-dimensional case is 
analogous to the [2 -dimensional] two-dimensional case. The 
10 specific method we have incorporated is a depth first search, 
presented by Aho, Hopecraft, and Ullman (A.V. Aho , [in " 1 J. E . 
Hopcroft, and J.D. Ullman. Design and Analysis of Computer 
Algorithms [, " section 5.2, published by 7 . Addison-Wes [t] ley, 
MA , 1974) . 

15 The result of partitioning at level M is the set of 



problems described as ( fP^ I i = 11 P,- J i=l , . . . , P-,- and j=l,...,N}, 
where N is the total number of hypothesis. The problems 
labeled { fP^li = 11 P< Ai=l , . . . ,p 2 } are the cases where 
rtheirl there is only one choice for the next observation at 
2 0 each scan and that observation could be used for no other 
track, i.e., it is a single isolated track. The hypothesis 
must be included in the solution set and no further processing 
is required . 
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As hypothesis are constructed the first element is used 
to refer to the track [id] ID. Any problem in the partitioned 
set which does not have shared costs in the first scan 
represent a case where a track could be extended in several 
5 ways but none of the ways share an observation with any other 
track. The required solution hypothesis for this case is the 
particular hypothesis with the maximum likelihood. For this 
case all likelihoods are determined as was described in 
Scoring and the maximum is selected. 

10 In addition to partitioning at the [M ] level M , 

partitioning is applied to each subsequent level M-l,...,2. 
For each problem that was not eliminated by [either by ] the 
prior selection^ partitioning is repeated, ignoring 
observations that are shared in the last set of observations. 

15 Partitioning recognizes that relaxation will eliminate the 
last constraint set and thus partitioning is feasible for the 
lower level problems that will result from relaxation. This 
process is repeated for all levels down to k=3 . The full set 
of partitionings can be performed in the Problem Formulation 

20 Step 310, prior to initiating the actual relaxation steps. 
The actual [2 -dimensional! two-dimensional solver used in step 
206 includes an equivalent process so no advantage would be 
gained by partitioning at the Jc=2 level. 
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There are two possible solution methods for the remaining 
problems. "Branch and Bound" as was previously described, or 
the relaxation method that this invention describes. If any 
of the partitions have 5-10 possible tracks and less than 50 
5 to 20 hypotheses, then the prior art "Branch and Bound" 
algorithm generally executes faster than does the relaxation 
due to its reduced level of startup overhead. The "Branch and 
Bound" algorithm is executed against all remaining M level 
problems that satisfy the size constraint. For the remainder 

10 the Relaxation algorithm is used. The scheduling done in 
Problem Formulation allows each Problem Formulation 
[ [3] il. [2] ] 51 cost matrix resulting from the first step of 
relaxation to be partitioned. The resulting partitions can be 
solved by any of the four methods [,] 4 isolated track direct 

15 inclusion, isolated track tree evaluation, small group "Branch 
and Bound" or an additional stage of relaxation as has been 
fully described. 

The partitioning after each level of Lagrangian 
Relaxation is effective because when a problem is relaxed, the 

2 0 constraint that requires each point to be assigned to only one 
track is eliminated (for one image at a time) . Therefore, two 
tracks previously linked by contending for one point will be 
unlinked by the relaxation which permits the point to be 
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assigned to both tracks. The fruitfulness of partitioning 
increases for lower and lower levels of relaxation. 

The following is a more detailed description of the 
partitioning method. Its application at all but the [M ] level 
5 M depends upon the relaxation process described in this 
invention. The recursive partitioning is therefore a distinct 
part of this invention. The advantage of this method is 
greatly enhanced by the sparse cost matrix resulting from 
tracking problems. However the sparse nature of the problem 

10 requires special storage and search techniques. 

A hypothesis is the combination of the cost element c n the 
selection variable z n and all observations that made up the 
potential track extension. It can be written as, 

h n ={c n ,z n , { [o nk =o k Jk=l]o n ^Jjc=l, . . . , [Mi] ^ e {1, . . . ,N k ) } } , i.e. , 

15 cost, selection variable and observations in the hypothesis. 
l"n£l Here n e {l,...,ivl, where N is the total number of 
hypotheses in the problem. While the cost and assignment 
matrices were previously referenced, these matrices are not 
effective storage mechanisms for tracking applications. 

20 Instead the list of all hypothesis and sets of lists for each 
dimension that reference the hypothesis set are stored. The 
hypothetical set in list form is: 

73 
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= For each dimension rk— 1 1 Jc=l , . . . j_M_ there exists a set of 
lists, with each list element being a pointer to a particular 
5 hypothesis: 

i=N kit k=M 



where N k is = § number of hypothesis containing the [i-th]i th 
observation from scan k. This structure is a multiply linked 
list in that any observation is associated with a set of 

10 pointer to all hypothesis it participates in, and any 
hypothesis has a set of pointers to all observations that 
formed it. (These pointers can be implemented as true 
pointers or indices depending upon the particular programming 
language utilized.) 

15 Given this description of the matrix storage technique 

then the partitioning technique is as follows: Mark the first 
list L k± and follow out all pointers in that list to the 
indicated hypothesis h Pk , for i=l, . . . , [N ki ] . Mark all located 
hypothesis, and for each follow pointers back the particular 

20 [L ok ] L k f° r ^=1/...^. Those L's if not previously marked get 
marked and also followed out to hypothesis elements and back 
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to L's. When all such L's or h x s being located are marked, 
then an isolated sub-problem has been identified. All marked 
elements can be removed from the lists and stored as a unique 
problem to be solved. The partitioning problem then continues 
5 by again starting at the first residual [L Jset^. When none 
remain, the original problem is completely partitioned. 

Isolated problems can result from one source track having 
multiple possible extensions or from a set of source tracks 
contending for some observations. Because one of the indices 
10 of k (in our implementation it is k=l) indicates the source 
track then it is possible to categorize the two problem types 
by observing if the isolated problem includes more than one 
L-list from the index level associated with tracks. 

15 1.4. Subsequent Track Prediction 



As noted above, after the points are assigned to the 
respective tracks, some action is taken such as controlling 
take-offs and landings to avoid collision, advising airborne 
20 aircraft to change course, warning airborne aircraft of an 
adjacent aircraft in a commercial environment, or aiming and 
firing an anti-aircraft (or anti-missile) missile, rocket or 
projectile, or taking evasive action in a military 
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environment. Also, the tracking can be used to position a 



applications, usually the tracks which have just been 
identified are extended or extrapolated to predict a 
5 subsequent position of the aircraft or other object. The 
extrapolation can be done in a variety of prior art ways; for 
example, straight line extensions of the most likely track 
extension hypothesis, parametric quadratic extension of the x 
versus time, y versus time, etc. functions that are the result 

10 of the filtering process described earlier, least square path 
fitting or Kalman filtering of just the selected hypothesis. 

The process of gating, filtering, and gate generation as 
they were previously described require that a curve be fit 
through observations such that fit of observations to the 

15 likely path can be scored and further so that the hypothetical 
target future location can be calculated for the next stage of 
gating. In the implementation existing, quadratic functions 
have been fit through the measurements that occur within the 
window. A set of quadratics is used, one per sensor 

2 0 measurement. To calculate intercept locations these 

quadratics can be converted directly into path functions. 
Intercept times are then calculated by prior methods based 
upon simultaneous solution of path equations. 



robot to work on an object. 



For some or all of these 
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The use of the fitted quadratics is not as precise as 
more conventional filters like the Kalman Filter. They are 
however much faster and sufficiently accurate for the relative 
scores required for the Assignment problem. When better 
5 location prediction is required then the assignment problem is 
executed to select the solution hypothesis and based upon the 
observations in those hypothesis the more extensive Kalman 
filter is executed. The result is tremendous computation 
savings when compared with the Kalman Filter being run on all 

10 hypothetical tracks. 

Based on the foregoing, apparatus and methods have been 
disclosed for tracking objects. However, numerous 

modifications and substitutions can be made without deviating 
from the scope of the present invention. For example, the 

15 foregoing 1 functions ¥ can also be defined as recursive 

approximation problem in which several values of higher order 
u k values are used to eliminate the higher than approximation 
characteristic of the [<$> k 1 function The hill climbing of 

the [W 1 function ¥ can be implemented by a high order hill 

20 climbing using the enhanced [<t> k 1 function <£l. Although the 
result would not be as efficient it seems likely that the 
method would converge. Therefore, the invention has been 
disclosed by way of example and not limitation, and reference 
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should be made to the following claims to determine the scope 
of the present invention. 

II . [ ] _AN ALTERNATIVE EMBODIMENT OF THE MULTI- 

5 DIMENSIONAL ASSIGNMENT SOLVING 

PROCESS 

The following description provides an alternative second 
embodiment of the multi -dimensional assignment solving process 
of section 1.1. However, in order to clearly describe the 
10 present alternative embodiment, further discussion is first 
presented regarding the data assignment problem of 

partitioning measurements into tracks and false alarms and, in 
particular, the representation of multi -dimensional assignment 
problems . 

15 

I I . 1 . Formulation of the Assignment Problem 
The goal of this section is to explain the formulation of 
2 0 the data association problems, and more particularly multi- 
dimensional assignment problems, that govern large classes of 
problems in centralized or hybrid centralized-sensor level 
multisensor/multitarget tracking. The presentation is brief; 
technical details are presented for both track initiation and 
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maintenance in ±A.&. Poore[, 1. Multi \ - 1 dimensional assignment 
formulation of data association problems arising from 
multitarget tracking and multisensor data fusion [, ]^ 
Computational Optimization and Applications, 3[ (1994), pp. 
5 ] : 27-57 , 1994 ) for non-maneuvering targets and in [Ingemar J. 
Cox, Pierre Hansen, and Beta Julesz, eds., Partitioning Data 
Sets, 1 (A.B. Poore . Multidimensional assignments and 
multitarget tracking: Partitioning data sets. In P. Hansen, 
I.J. Cox, and B. Julesz, editors, DIMACS Series in Discrete 

10 Mathematics and Theoretical Computer Science , volume 19, pages 
169-198, Providence, R.I., 1995 . American Mathematical 
Society[, Providence, R.I., v. 19, Feb. 1995, pp. 169-198]! 
for maneuvering targets. Note that the present formulation 
can also be modified to include target features (e.g., size 

15 and type) into the scoring process 154. 

The data assignment problems for multisensor and 
multitarget tracking considered here are generally posed as 
that of maximizing the posterior probability of the 
surveillance region (given the data) according to 



20 



[[5.1] 



Maximize 



i P(r=y\Z M ) 



(2.1) 



p(r=y°\z M ) 



where Z M represents M data sets, y is a partition of 
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indices of the data (and thus induces a partition of the 

data) , T* is the finite collection of all such partitions, r 

is a discrete random element defined on r*, y° is a reference 

partition, and P( [r=vlz M 1 F=v\z M ) is the posterior probability 

5 of a partition y being true given the data Z". The term 

partition is defined below. Consider M observation sets 

r k , N. ] f if iWjt 

Z{k) (k=l, . . . ,N) each of N k observations [{ z ^ I . K - \z\ ) , and 

x k i k -l I k h k =i 

let Z M denote the cumulative data set = defined by 
[ and Z M = {Z(l) , . . . ,Z(M) }, 

10 

[5.2] 



Z(k) ={z- k J k i and Z m = {z(l),...,Z(M)}, 



(2.2; 



] respectively . [ ] In multisensor data fusion and multitarget 
tracking the data sets Z(k) may represent different classes of 

15 objects, and each data set can arise from different sensors. 
For track initiation the objects are measurements that must be 
partitioned into tracks and false alarms. In the formulation 
of track extensions a moving window over time of observations 
sets is used. The observation sets will be measurements which 

20 are: (a) assigned to existing tracks, (b) designated as false 
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measurements, or (c) used for initiating tracks. However, 
note that alternative data objects instead of observation sets 
may also be fused such as in sensor level tracking wherein 
each sensor forms tracks from its own measurements and then 
5 the tracks from the sensors are fused in a central location. 
Note that, as one skilled in the art will appreciate, both 
embodiments of the present invention may be used for this type 
of data fusion. 

The data assignment problem considered presently is 

10 represented as a case of set partitioning defined in the 
following way. First, for notational convenience in 

representing tracks, a zero index is added to each of the 
index sets in [ [5] X2. 2 [] ] J =/ a gap filler z k 0 is added to each of 
the data sets Z(k) in [ [5] Ag. 2 [] ] 1, and a "track of data" is 

15 defined as = ^ where i k [and zj can now assume the values of 
0 [ and z, respectively] . A partition of the data will refer 
to a collection of tracks of data or track extensions wherein 
each observation occurs exactly once in one of the tracks of 
data and such that all data is used up. Note that the 

20 occurrence of the gap filler is unrestricted. The gap filler 
Zq serves several purposes in the representation of missing 
data, false observations, initiating tracks, and terminating 
tracks. Note that a "reference partition" may be defined 
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which is a partition in which all observations are declared to 
be false. 

Next under appropriate independence assumptions it can be 
shown that : 
5 [[5.3]] 



P(r- Y |Z") = £ = n [ L L ( 2.3) 



where L* is a likelihood ratio containing probabilities for 
detection, maneuvers, and termination as well as probability 
10 density functions for measurement errors, track initiation and 
termination. Then if c? I » M -^,^= -ln(L^ , Jf ) , 
[ [5.4] 

[ c. c, 1 . (2.4) 



[]-ln 

Us i ng 



P( [y]r=y\z M ) 



(i x , . . .,i M )6Y 



_P( [y ]T=Y |2 
[5.3] 

15 and the zero-one variable z x> .. iN =l if (ii,...,^) 6 y and 0] 

Using (2.3) and the zero-one variable 

z* 2 < jr =l. i£ (i ^ i^) £ Y' and °' otherwise, one can then 

write the problem [ [5] Jj.. 1 [] ] J. as the following M-dimensional 
assignment problem (as also presented in section Problem 

20 Formulation [[3]ll. [l]]4l with k=M) : 
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N. N, 



(a) [ ] Minimize 



(b) [ ] Subject to. 



i 1= o i M =o 



M, N. 



[]52---5J z i i =l / [ ] ii = i, . . . ,^1, 



i 2 =0 i M =o 



(C) [ ] [] E • • • E E • • • E *i i =lr (2.5) 

i 1 = 0 1^=0 i w n0 i M =0 1 " 

[ ] for = 1, . . . ,N k and k=2, . . . ,M-1, 
5 (d) [ ] [] E • • • E *i i =1, [ 1 i M = 1/ • • • 

(e) [ ] n z i r .,i/ {0 ' l} for all i^...,^ 

[ [5.5] 
] where c 0 ... 0 is arbitrarily defined to be zero. Here, each 
10 group of sums in the constraints \ \S] (2 . 5 M 1 ) (b) - 
[ [5] . 5 [] ] 1 (e) represents the fact that each non-gap filler 
observation occurs exactly once in a "track of data." One can 
modify this formulation to include multi^assignments of one, 
some, or all the actual observations. The assignment problem 
15 [ [5] 12.. 5 [] ] 1 is changed accordingly. For example, if z\ k is to 
be assigned no more than, exactly, or no less than n k k times, 
then the "=l" [in the appropriate one of the constraints 
[5 . 5] (b) - [5 . 5] (d) ] is changed to "£, =, £ n k ik , " respectively. 
Expressions for the likelihood ratios L k 1 /_ ±M [ such as 
2 0 [5.5] ] appearing in the costs c L i =-ln(L i 1 ) are well [ ] ^ 



83 



NUMERI .01USC2 

known. In particular, discussions of these ratios may be 
found in [ : 1 (A.B. Poore[, 1 . Multi [-] dimensional assignment 
formulation of data association problems arising from 
multitarget tracking and multisensor data fusion, 
5 Computational Optimization and Applications [ , 3 (1994), pp. 
27-57; and Ingemar J. Cox, Pierre Hansen, and Beta Jules z, 
eds w Partitioning Data Sets, 7 , 3:27-57, 1994; A.B. Poore . 
Multidimensional assignments and multitarget tracking: 
Partitioning data sets. In P. Hansen, I.J. Cox, B. Julesz, 

10 editors , DIMACS Series in Discrete Mathematics and Theoretical 
Computer Scienc e, volume 19, pages 169-198, Providence, R.I., 
1995 . American Mathematical Society [, Providence, Rll. [I., v. 
19, 1995, pp. 169-198. ] Furthermore, the likelihood ratios 
are easily modified to include target features and to account 

15 for different sensor types. Also note that in track 
initiation, the M observation sets provide observations from 
M sensors, possibly all the same. Additionally note that for 
track maintenance, a sliding window of M observation sets and 
one data set containing established tracks may be used. In 

20 this latter case, the formulation is the same as above except 
that the dimension of the assignment problem is now M+l. 



I I . 2 . Overview of the New Lagrangian 
Relaxation Algorithms 
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Having formulated an M-dimensional assignment problem 
[ [5] Jj^ 5 [] ] 1, we now turn to a description of the Lagrangian 
relaxation algorithm. Subsection II. 2.1 below presents many 
of the relaxation properties associated with the relaxation of 
5 an n -dimensional assignment problem to an jn-dimensional one 
via a Lagrangian relaxation of n-m sets of constraints wherein 
[m<1 m < n < M and preferably in the present invention 
embodiment 

[n-m>ll n-m >1 . Although any n-m constraint sets can be 
10 relaxed, the description here is based on relaxing the last n- 
m sets of constraints and keeping the first m sets. Given 
either an optimal or suboptimal solution of the relaxed 
problem, a technique for recovering a feasible solution of the 
n-dimensional problem is presented in subsection II. 2. 2 below 
15 and an overview of the Lagrangian relaxation algorithm is 
given in subsection II. 2. 3. 

The following notation will be used throughout the 
remainder of the work. Let M be an integer such that [M^3]£f 
£3 and let ne{3,...,M}. The n-dimensional assignment problem 
20 is 
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(a) Minimize: v n (z) =[ J^z2 • • • zJ c ± l ...i n z i 1 . ,.i n 



i x =0 i„=0 



(b) [ 
(c) 

(d) 
(e) 



] Subject [to: 



ij = 1, . . . 

for i*. = 1, . . . ,N k and k=2, . . . ,n-l, 
i n = I/..., N n , 

for all i lf . . . , i n ] 



10 



15 (el 



[6.1] ]To 



2^/ ••'2^/ Z i 1 ...i n = ^-r 1 i = 1' • • • ' ®\i 



^Lj ••• £ S z i,...i 1 

for 1^ = 1, . . . , A/^ and k = 2, . . . , n-1, 
s i 1 ...i n e {°' 1 } for all i 1 ,...,i n . 



(2.6) 



To ensure that a feasible solution of r M (2 . 6 T . 11 1 ) always 
exists, all variables with exactly one nonzero index (i.e., 
20 variables of 

the form z n 0 . . . oik0 . . . 0 for 2**0) are assumed free to be assigned 
and the corresponding cost coefficients are well-defined. 
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[ ] a _II.2.1. The Lagranaian Relaxed Assignment 

Problem 

The n-dimensional assignment problem [ M (2 . 6 [ . 1] 1 ) has n 
sets of constraints. A {N k +1) -dimensional multiplier vector 
5 associated with the [k-th]^ constraint set will be denoted by 
u k = (u£, u k lf . . . , u^ k ) T with u£=0 and k=l, . . . ,n. The n-dimensional 
assignment problem \ M (2 . 6 [ , 1] ] j_ is relaxed to an m- 
dimensional assignment problem by incorporating n-m of the n 
sets of constraints into the objective function 
10 [ n (2 . 6 r . 11 1 ) (a) . Although any constraint set can be relaxed, 
the description of the relaxation procedure for [ [] (2 . 6 [ . 11 1 ) 
will be based on the relaxation of the last n-m sets of 
constraints. The relaxed problem is 
[^(u^ 1 , . . . ,u n ) = Minimize 0^ (z n ;u m+1 , . . . , u n ) 

15 = 
+ 

Subject to: i x = 1,...,N 1# 

2 0 for i k = 1 , . . . , N k and k=2 , . . . , m, 

for all i 1# . . . , i n . ] 
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mn * ' 



,u n )= Minimize Q mn (z n ; u m+1 , . . . , u") 

S " ' Ci l--- i n Zi l--- i n 



(2.7) 



i^o 



£ £ u*. 



Z> ••• L~i £-J ' ' ' Z~i z i....i„ 



[6.2] ] 



= c i...i n + £ u i z i...i n ~ £ £ u i* 

i.=0 1=0 1 " Jt-a+1 * 1 n 1 * 

1 n L J 

y ; • • • / ^ ^ i lt . . i n ~ 1 ' 1 ~ ^" ' * " * r ^1 ' 



Subject To 



5-/ • • * • • • 5-) z ±,...± n i 

for 1^ = 1, . . . , N k and k = 2, . . . ,m, 
^...i^fO'l} for all i 1# ...,i n , 



wherein we have added the multiplier u£=0 for notational 
convenience. Thus, the multiplier with u$=0 is fixed. 

An optimal (or suboptimal) solution of 
[ [6] J2. . [2] ] 21 can be constructed from that of an m- 
dimensional assignment problem. To show this, define 
for each an index ( j m+1/ . . . , j n ) = (j m +i([ii 
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(] i 2 , . . . , i J -L (i w . . . , iJ ) and a 

new cost function /"c^.^J by 

[ (j m+ w • • • / j n )=arg min{c 1 ... in +u lt |i k =0 f . . . ,N k and 

k=m+l , . . . , n} 



m+1' 



irgmin \c[ ' ->#1 + 2^ u£ | 1^ = 0, . . . , and /c = in+l, . . . , nL (2. 



k...o £ < 

[ 1 " Jt=a+1 k 

i + £ "A. for (iw...,iJM0,...0), 



[6.3]] 

[If (j m+1 , . . . , j n ] 

C*. . .o = 2 • • • £ Minimum { 0, c 0 " . >01 . . j + £ u/ 1. 



If ( ~L , . . . , 7 J is not unique, choose the smallest such 
amongst those (n-m) -tuples with the same j m+1 choose 
the smallest j m+2 , etc., so that (j BH . 1 ,-..,j n ) is uniquely 
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defined.) Using the cost coefficients defined in this 
way, the following m- dimensional assignment problem is 
obtained: 

g^iif* 1 , . . . ,u n ) = Minimize ^(z"; u» ... , u n ) = v m (z m ) 
• • • 1-, . . . i m ,.i m 

- . 1 ml m 

{ ] Subject to: Y, . . . z" . . #i =1 , [ ] i 2 = 1, . . . ,N lf 

• n 'A 1 " " " O 



[ ] 



5-^ " • • ^L/ " " " 5-) Z i 1 ...i_ ^' 



[ ] for i k = 1, . . . ,1^ and Jc=2, , . . 



5-^ • • • > %i . . . i If — / • • • / 

[ ] <...i e fo f l} for all i lf .. .,i„. 

[ [6.4] 
As an aside, observe that any feasible solution z n of J2^6 . [1]^- 
yields a feasible solution z m of [6U£. ] 4] 9^ via the 
construction 



Thus, the jn-dimensional assignment problem [ [6] X2. [4] ] '9X has 
at least as many feasible solutions of the constraints as the 
original problem \ \] (2 . 6 [ . 1] ] ) . 
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The following Fact A. 2 has been shown to be true. It 
states that an optimal solution of Problem Formulation 
[[6]J2.. [2]]XL can be computed from that of Problem Formulation 
[ [6] [4] ] £1. Furthermore, the converse to this fact is 
5 provided in Fact A. 3. Moreover, if the solution of either of 
these two problems [ [6]ig. [2] ] 71 or [ [6]ig. [4] ] 91 is 6-optimal 
(i.e., the objective function associated with the suboptimal 
solution is within " E" of the objective function associated 
with the optimal solution) , then so is the other. 
10 Fact A. 2. Let vJ" be a feasible solution to problem 

[ [6]J£. [4] ] 91 and define w" by 

Wi 1 ...i„=wT J ...i„ if • • -i n ) = • • • J n ) and 

(i lf . . . ,ij * (0, . . . , 0) 

15 w? 2 ...i n =0 if d m+1 , . . .i n ) * (J m+ i/ • • • 3n) and (2.10) 

(i 2/ . . . , ij * (0, . . . , 0) 

<.. ow .. in =l if cS... 0W ... Ja + £0 
<..ow...i n =0 if cs... 0im+1 ... in + >0^ 

[[6.5]] 

20 Then w" is a feasible solution of the Lagrangian relaxed 
problem [ [6] J£. [2] ] 71 and 

<P mn (w n ;u m + 1 ,...,u n ) 



91 
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If, in addition, vf is optimal for [[6]A2. [4]] 91, then w" is an 
optimal solution of [ [6] 12 . [2] ] 71 and 

^„(u" tJ ,...,u") 

" •» (um+1 u n > "E [ £=m + l£ VO u ^ J t < . 



5 With the exception of one equality being converted to an 

inequality, the following Fact is a converse of Fact A. 2 and 
has also been shown to be true. 

Fact A. 3. Let w" be a feasible solution to problem 
[ [6]12. [2] ] 71 and define vf by 

10 wf i = [] Jl ^...in for • • -/iJMO, . . .,0)^ 

wS... 0 =0 if (i 2 , . . .,i m ) = (0, . . .,0) and 

Co...ow..i n + >0 = for a22 iJ • (2.11) 

<... 0 =1 if (ii/ • • • / iJ = (0/ • ■ • / 0) an & 
&o...oi m¥l ...i n + ^° for some (i m+a / • . wi n ) • 
15 [[6.6]] 

Then w™ is a feasible solution of the problem 

[ [61^2,. [4] ] 21 and , i ^(w";u ,tl ,... ( u") * 
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If, in addition, w" is optimal for [[6]J2:. [2]]2l/ then W" is an 
optimal solution of [ [6] 12 . [4] ] 11, [0^ (w n ;u m+1 , . . . ,u n ) = , and 
<!>„, (u m+1 ,u n ) = . ] 



[ 



<P mn (w";u m +\... f u") =<f mn {w m ;u m + 1 ,...,u n ) - £ u k i f and 

*iu» + \...,u n ) =£ B <u- i ,...,u n ) - £ £ u}. 



k=m+l i k =0 k 



5 ] II. 2. 2. The Recovery Procedure 

The next objective is to explain a recovery procedure or 
method for recovering a solution to the n-dimensional problem 
of r n (2 . 6 f . Ill ) from a relaxed problem having potentially 
substantially fewer dimensions than r M (2 . 6 [ . 1] ]j. Note that 

10 this aspect of the alternative embodiment of the present 
invention is substantially different from the method disclosed 
in the first embodiment of the multi [-] dimensional assignment 
solving process of section 1.1 of this specification. 
Further, this alternative embodiment provides substantial 

15 benefits in terms of computational efficiency and accuracy 
over the first embodiment, as will be discussed. Thus, given 
a feasible (optimal or suboptimal) solution w™ of [[6]JJ2. [4] ]9_L 
(or vf 1 if Problem Formulation [[6)iZ. [2]] 7± is constructed via 
Fact A. 2), an explanation is provided here regarding how to 
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generate a feasible solution z n of [ [] (2 . 6 [ . 1] ] j_ which is close 
to w™ in a sense to be specified. We first assume that no 
variables in [ [] (2 . 6 [ . 1] ] ) are preassigned to zero (this 
assumption will be removed shortly) . The difficulty with the 
5 solution vf 1 is that it need not satisfy the last n-m sets of 
constraints in I M (2 . 6 [ . 1] ] ) . Note, however, that if is an 
optimal solution for [ [6] [4] ] 9X and (constructed as in 
Fact A. 2) satisfies the relaxed constraints, then w" is optimal 
for I" M (2 ■ 6 [ ■ 1] ] ) . The recovery procedure described here is 

10 designed to preserve the 0-1 character of the solution of 
[[6]I2. [4] ]£1 as far as possible. That is, if w? 2 _ in =l and i^O 
for at least one 1=1,..., m, then the corresponding feasible 
solution z n of f M (2 , 6 [ . 1] ]1 is constructed so that z2 1<i<in =l for 
some (i m+1 , . . . , i n ) . Note that by this reasoning, variables of 

15 the form z?. . .oi m+I . ..i„ can be assigned to a value of 1 in the 
recovery problem only if wg_ 0 =l. However, variables .oi m+1 .. .± n 
will be treated differently in the recovery procedure in that 
they can be assigned 0 or 1 independent of the value wJJ... 0 . 
This increases the feasible set of the recovery problem, 

20 leading to a potentially better solution. 

Let [{(ii), . . . , i^)} N ° , ]| (i/, ...fifl) [ W ° be an enumeration of 
l mj=ii x ) j-i 

indices of w m (or the first m indices of w" constructed as in 
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m 

Fact A. 2) such that [w. j . j = l]w m j , } = 1 and ( i{, . . . , i J m ) * 
(0,...,0). Set (i° lf . . .,i°) = (0, . . .,0) for j = 0 and define 
[ for i k =0,...,N k ; k=m+l,...,n; j=0,...,N 0 . 

[6.7] 

n-m+1 _ ^ n 

Let Y] C ji m+1 ...i n " <V ± 3 ± i = ; 



(2.12) 

for i „=0 ..... N i . ; k=m+ 1 , . . . , n ; i =0 . . . . . N n . 



Let y denote the solution of the (n-m+1) -dimensional 
assignment problem: 

Minimize 



Subject to 
15 (2 . 13) 



N„ W-., Nu.i Ni.^ N, 



□ 2, Is ••• 2^ 2^ •••2^y J1 ^..i^i/ 

j=0 i„ n = 0 1^=0 i n =0 -» - 

for i k =l, . . . ,N k and k=m+l, . . . ,n-l, 

[]E E ••■ E y« i =1 ' i ] i„=i, . . . ,w ra , 

[ly JW"*n 6 {0 ' 1} f ° r a21 ^ i "' 

[[6.8]] 

The recovered feasible solution z n of f f 1 (2 . 6 f . 11 1 ) 
corresponding to the multiplier set {u m , . . .,u n } is then defined 
by 
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n 



1, if . . . , ij = (il, . . . , i£) for some 

j=0, . . . ,N 0 and Y . ..i=l, (2.14) 

0, otherwise. 



[ [6.9] 

]This recovery procedure is valid as long as all cost 

coefficients c" are defined and all zero-one variables in z n 

are free to be assigned. Note that modifications are 

5 necessary for sparse problems. If the cost coefficient 
n 

[c . j . j . . ] c n j j [ ] is well defined and the zero- 

1 1' ' 'Vm+l' • nn 2 i 
one variable [z.j .j. . ] c n j j [ ] is free to be 

1 1* ' 'Vm+l' # m:L n 1 i-- 1 * 1 *+i'- 1 n 

assigned to zero or one, then define 

n-m+1 n , n . m+1 n . 

j = c.j .j. . ] c-n i = c j i [ ] as m 

10 [[6]^. r7l 1 12) with being free to be assigned. Otherwise, 

is preassigned to zero. To ensure that a feasible solution 
exists, we now only need ensure that the variables are free 
for j=0, 1, . . . f N 0 . (Recall that all variables of the form 
z o...i k ...o are free (to be assigned) and all coefficients of the 

15 form cS_ ti 0 are well defined for k=l , . . . , n. ) This is 

accomplished as follows. If the cost 

n 

c[c.j.j -in r\]c n i 1 i oefficient is well defined and 

[z.j.j . j n r\]z n i i i is free, then define 

1 z m n 

r n-m+1 . . . , n-m+1 n . , . 

[C jO...O = C ±iq...il Q...0 ]c JO...o = c iih i. ..^o .. .o Wlth = bein 9 
20 free. Otherwise, since all variables of the form z^. . . .o are 
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known feasible and have well-defined costs, put 

r n-m+l _ ^ n . . n-m+l _ ^ n 

Lc jO...O " A( C 0...0ij0...0 JC J°---° ^ C o...oi^o...o * 



II. 3. The Multi -Dimensional Lagrangian Relaxation 

5 Algorithm for the n-Dimensional Assignment 

Problem 

Starting with the M-dimensional assignment problem 
r n (2 . 6 [ . 1] ] Jl, i.e., n=M, the algorithm described below is 
recursive in that the M-dimensional assignment problem is 

10 relaxed to an m-dimensional problem by incorporating (n-m) 
sets of constrains into the objective function using 
Lagrangian relaxation of this set. This problem is maximized 
with respect to the Lagrange multipliers, and a good 
suboptimal solution to the original problem is recovered using 

15 an (n-in+1) -dimensional assignment problem. Each of these two 
(the /n-dimensional and the (n-m+l) -dimensional assignment 
problems) can be solved in a similar manner. 

More precisely, reference is made to Figure 8 which 
presents a flowchart of a procedure embodying the multi- 

20 dimensional relaxation algorithm, referred to immediately 
above. This procedure, denoted MUL T I _D I M_RE LAX in Figure 8, 
has three formal parameters, n, PROB_FORMULAT I ON and 
ASSIGNMT_SOLUTION, which are used to transfer information 
between recursive instantiations of the procedure. In 
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particular, the parameter, n t is an input parameter supplying 
the dimension for the mult i -dimensional assignment problem (as 
in [ [1 (2 . 6 [ . 1] ] ) ) which is to be solved (i.e., to obtain an 
optimal or near-optimal solution) . The parameter, 

5 PROB_FORMULAT I ON , is also an input parameter supplying the 
data structure (s) used to represent the n-dimensional 
assignment problem to be solved. Note that PROB_FORMULAT I ON 
at least provides access to the cost matrix, [c n ] , and the 
observation sets whose observations are to be assigned. 

10 Additionally, the parameter, ASSIGNMT_SOLUTION, is used as an 
output parameter for returning a solution of a lower 
dimensioned assignment problem to an instantiation of 
MULT I _D I M__RELAX which is processing a higher dimensioned 
assignment problem. 

15 A description of Figure 8 follows. Assuming 

MUL T I _D I M_RE LAX is initially invoked with M as the value for 
the parameter n and the PROB_FORMULAT I ON having a data 
structure (s) representing an M-dimensional assignment problem 
as in [ [5] 12. 5 [] ] 1, in step 500 an integer m, r2^m<n1 2£m<n is 

20 chosen such that the constraint sets fori corresponding to 
dimensions m+l,...,n are to be relaxed so that an Tri- 
dimensional problem formulation as in [ [6] S2 . [2] ] 71 results. 
In step 504 an initial approximation is provided for 
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{uq +1 , . . . , Uq) . Subsequently, in step 508 the above initial 
values for { u™* 1 , . . . , u"} are used in an iterative process in 
determining {u 1 }IU+1 , . . . /= y 7 ]n ) which maximizes {^(u" 111 , . . . , u n ) 
where ru k eR Mk+1 1 u*eR M * +I and k=m+l , . . . , n} for a feasible solution 
5 w" subject to the constraints of Problem Formulation 
[ [6] 12. [2] ]£L- Note that by maximizing cP^ (u m+1 , . . . , u n ) ., some 
of the constraints that were relaxed are being forced to be 
satisfied and in so doing information is built into a solution 
of this function for solving the input assignment problem in 

10 " PROB__FORMULATION . " Further note that a non-smooth 

optimization technique is used here and that a preferred 
method of determining a maximum in step 508 is the bundle- 
trust method of Schramm and Zowe [as described in ]±H. Schramm 
and J. Zowe[, 1 . A version of the bundle idea for minimizing 

15 a non- smooth function: Conceptual idea, convergence analysis, 
numerical results [, ] . SI AM Journal on Optimization, 2[ 
(1992)], [pp]No. 1:121-152 , February, 1992) . This method, 
along with various other methods for determining the maximum 
in step 508, are discussed below. 

20 Also note that for [vc\>2] m >2 , a solution to the 

optimization problem of step 508 is NP-hard and therefore 
cannot be solved optimally. That is, there is no known 
computationally tractable method for guaranteeing an optimal 
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solution. Thus, there are two possibilities: either (a) allow 
m to be greater than 2 and use auxiliary functions similar to 
those disclosed in the first embodiment of the Jc-dimensional 
assignment solver 3 00 in section I to compute a near-optimal 
5 solution, or (b) make m-2 wherein algorithms such as the 
f orward/reverse auction algorithm of D.P. Bertsekas and D.A. 
Castanon [in their paper, ] (D.P. Bertsekas and D.A. Castanon. 
A forward/reverse auction algorithm for asymmetric assignment 
problems [, 1 . Computational Optimization and Applications, 1[ 

10 (1992), pp.]^ 277-298 , 1992) provides an optimal solution. 

If option (a) immediately above is chosen, then the 
auxiliary functions are used as approximation functions for 
obtaining at least a near-optimal solution to the optimization 
problem of step 508. Note that the auxiliary functions used 

15 depend on the value of m. Thus, auxiliary functions used when 
m=3 will likely be different from those for m=4 . But in any 
case, the optimization procedure is guided by using the merit 
function, <t> 2n ( u 3 , . . . , u n ) , which can be computed exactly via a 
[2-dimensional] two-dimensional assignment problem for guiding 

2 0 the maximization process. 

Alternatively, if option (b) above is chosen, then two 
important advantages result [ : ] , namely , the optimization 
problem of step 508 can be always solved optimally and without 
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using auxiliary approximation functions. 



Thus, better 



solutions to the original M-dimensional assignment problem are 
likely since there is no guarantee that when the non- smooth 
optimization techniques are applied to the auxiliary functions 
5 the techniques will yield an optimal solution to step 508. 
Furthermore, it is important to note that without auxiliary 
functions, the processing in step 508 is both conceptually 
easier to follow and more efficient. 



10 (u 1 , . . . , lT n ) [and w n ] is used in determining an optimal 
solution[, ] == w m to Problem Formulation [ [6]^2. [4] ]£L as 
generated according to [ [6] J^. [3] ] JLL and Fact A. 3. 

In step 516, the solution is used in defining the cost 
matrix c n " m+1 as in r f6l (2 . [7] ] 12) . Subsequently, if n-m+l=2, 

15 then the assignment problem r T61 (2 . T81 1 13) may be solved 
straightforwardly using known techniques such as forward/ 
reverse auction algorithms. Following this, in step 528, the 
solution to the f 2 -dimensional! two-dimensional assignment 
problem is assigned to the variable ASSIGNMT_SOLUTION and in 

20 step 532 ASSIGNMT_SOLUTION is returned to a dimension three 
level recursion of MULT I _D I M_RELAX for solving a [3- 
dimensional] three-dimensional assignment problem. 



Subsequently, in step 512 of Fig. 8, the solution 
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Alternatively, if in step 520, n-/n+l>2, then in step 536 
the data structure (s) representing a problem formulation as in 
[ [6] (2 . r 8 1 1 13 ) is generated and assigned to the parameter, 
PROB_FORMULATION. Subsequently, in step 54 0 a recursive copy 
5 of MULTI_DIM_RELAX is invoked to solve the lower dimensional 
assignment problem represented by PROB__FORMULATION. Upon the 
completion of step 540, the parameter, ASSIGNMT_SOLUTION, 
contains the solution to the [n-m+1 dimensional! (n-ra+1) - 
dimensional assignment problem. Thus, in step 544, the [n- 

10 m+ll (n-m+1) th solution is used to solve the n-dimensional 
assignment problem as discussed regarding equations 
[[6112. r 9 1 1 14) . Finally, in steps 548 and 552 the solution to 
the n-dimensional assignment problem is returned to the 
calling program so that, for example, it may be used in taking 

15 one or more actions such as (a) sending a warning to aircraft 
or sea facility; (b) controlling air traffic; (c) controlling 
anti-aircraft or anti-missile equipment; (d) taking evasive 
action; or (e) surveilling an object. 

20 II. 3.1. Comments on the Various Procedures 

Provided by Figure 8 

There are many procedures described by Figure 8 . One 

such procedure is the first embodiment of the multi [- 

25 ] dimensional assignment solving process of section 1.1. That 
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is, by specifying m=n-l in step 500, a single set of 



constraints is incorporated [ is] into the objective function 
via the Lagrangian problem formulation, = resulting in an [m=n-l 
5 dimensional! (m=n-l) -dimensional problem. The relaxed problem 
is subsequently maximized in step 512 with respect to the 
corresponding Lagrange multipliers and then a feasible 
solution is reconstructed for the n -dimensional problem using 
a two-dimensional assignment problem. The second procedure 

10 provided by Figure 8 is a novel approach which is not 
suggested by the first embodiment of section 1.1. In fact, 
the second procedure is somewhat of a mirror image of the 
first embodiment in that n-2 sets of constraints are 
simultaneously relaxed, yielding immediately an [m=2 

15 dimensional! (m=2 ) -dimensional problem in step 512. Thus, a 
feasible solution to the n-dimensional problem is then 
recovered using a recursively obtained solution to an [n-1 
dimensional! (n-1 ) -dimensional problem via step 540. In this 
case, the function values and subgradients of & 2n (u 3 , . . . , u 11 ) of 

20 step 508 can be computed optimally via a two [ ] ^dimensional 
assignment problem. The significant advantage here is that 
there is no need for the merit or auxiliary functions, W kt as 
required in the first embodiment of the multi [-] dimensional 



constraints is relaxed in step 508. 



Thus, one set of 
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assignment solving process of section I.l above and also there 
is no need for the more general merit or auxiliary functions 
Wxn, as discussed in subsection II. 4. 2 below. Further, note 
that all function values and subgradients used in the non[- 
5 ] smooth maximization process are computed exactly (i.e., 
optimally) in this second procedure. Moreover, problem 
decomposition is now carried out for the n-dimensional 
problem; however, decomposition of the [n-1 dimensional! (n- 
1) -dimensional recovery problem (and all lower recovery 

10 problems) is performed only after the problem is formulated. 

Between these two procedures are a host of different 
relaxation schemes based on relaxing n - m sets of constraints 
to an /n-dimensional problem ( f2<m<n1 2<m<n ) , but these all have 
the same difficulties as the procedure for the first 

15 embodiment of section I.l in that the relaxed problem is an 
NP-hard problem. To resolve this difficulty, we use an 
auxiliary or merit function as described in subsection 

II. 4. 2 below. (The notation W k was used for W n . ltTl in section 
I.l.) For the case m<n-l, the recovery procedure is still 

2 0 based on an NP-hard 

( rn-m+l-dimensionall n-m+1) -dimensional assignment problem. 
Note that the partitioning techniques similar to those 
discussed in Section 1.3 may be used to identify the 
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assignment problem with a layered graph and then to find the 
disjoint components of this graph. In general, all relaxed 
problems can be decomposed prior to any non[-] smooth 
computations because their structure stays fixed throughout 
5 the algorithm of Fig. 8. However, all recovery problems 
cannot be decomposed until they are formulated, as their 
structure changes as the solutions to the relaxed problems 
change. [ ] 

10 

II. 4. Details and Refinements Relating to the 

Flowchart of Figure 8 

Further explanation is provided here on how various steps 

15 of Figure 8 are solved. Note that the refinements presented 

here can significantly increase the speed of the relaxation 

procedure of step 508. 

II. 4.1 Maximization of the Non-smooth Function 
20 A ndi"* 1 ,--.,"") 

One of the key steps in the Lagrangian relaxation 
algorithm in section II. 3 is the solution of the problem 
[ ] 
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= Maximize {^(u^ 1 , . . . f u n ) : ru k eR Mk+1 1 u k el M * +I ; [k=m+l]k= 
/n+1 , . . . , n} [ 

[7]* J2.1 n ] 51 

5 where u£=0 for all k=m+l , . . . , n. The evaluation of 

^ (u™* 1 , • . . / u n ) requires the optimal solution of the 
corresponding minimization problem [[6]12. [2] Z n varies for 
each instance of (u m+1 , . . . , u n ] 7) . [ ] The following discussion 
provides some properties of these functions. 

10 Fact A. 4. Let u m+1 , . . . , u n be multiplier vectors associated 

with the (m+2) st through the n th set of constraints 
[ [] (2.6 Ml ] ) , let <t> m be as defined in [ [6] 12. [2] ] 21, let 
V n {zJ ]n ) be the objective function value of the n-dimensional 
assignment problem in equation [ M (2 . 6 [ . 1] 3 ) , let z n be any 

15 feasible solution of r M (2 . 6 [ . 1] ] 1, and let i 7 ]n be an optimal 
solution of I" M (2 . 6 \ . 11 1 ) . Then, {u m+1 , . . . , u n ) is piecewise 
affine, concave and continuous in { u m+1 , . . . , u"} and 

0~- (u m+J u n ) £ v n (F n ) z v n (z n L^ = e _ s _ 

20 (2 . 16 ) 

[[7.2a]] 

Furthermore, 
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0 m . ln (u m , U m+1 , ...,U n ) <L 0^ (U m+1 , U nI > 

[7.2a] 

for all ukeRMk+1 with u=0 and k=m, . . . , n . ] j (2 17) 

for all u*6R M * +1 with u£=0 and k=m n. 



5 



Most of the procedures for non[-] smooth optimization are 
based on generalized gradients called subgradients , given by 
the following definition. 

Definition . At u= {u m * x , . . . , u n ) the set [ 5<P„„] d<P m „(u) is called 
10 a subdifferential of 0^ and[ is] defined by 
[SOU 

§0 mn (u)={ rael Mm+1+1 xl a € R*^ +I x . . . rxR Mn+1 l$„J xR** +I l <p_ ( w ) [-tj- 

grnn(u) [<q T ]^( [ W- U ) VwER Mm+1 + 1 X . . . xR Mn+1 } 

[7.3] 

15 A vector qeSO^di) 1 w-u) (2.18) 

Vw E R M ^ +i x. . .xR^ +1 ) . 

A vector a 6 d& m „(u) is called a subgradient . 

If z n is an optimal solution of [ [6] J^. [2] ] 2X computed 
20 during evaluation of ^(u) , differentiating with respect to 
u i n yields the following [i n -th]i^ component of a subgradient 
9 of <?U(u) 

k n 
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i^O i w -0 i^O i n =0 

for i k =l, . . . and =771+2, ... ,77. 

[[7.4]] 

If z n is the unique optimal solution of r [ 6 1 ( 2 . [ 2 ] ] 7 ) , 
5 [50^ (u) ={g} , and 3^] d^ Ju) = { a] and^ is diff erentiable at u. 
If the optimal solution is not unique, then there are finitely 
many such solutions, say z n (1) , . . . , z n (K) . Given the 

corresponding subgradients , gr 1 /.../^*/ the subdif f erential 
[5<Z>]d&(u) is the convex hull of {gr 1 , . . . , gr*} . 

10 

II. 4. 2 Mathematical Description of 

a Merit or Auxiliary Function 

For real-time needs, one must address the fact that the 

non-smooth optimization problem of step 508 (Fig. 8) requires 

15 the solution of an NP-hard problem for m>2 . One approach to 

this problem is to use the following merit or auxiliary 

function to decide whether a function value has increased or 

decreased sufficiently in the line search or trust region 

methods : 
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r • • • 9 




u n ) = 



(2.20) 




jn+l 



i • • • r 



if jn = 2 

or ([6] 2. [2] 7) is solved optimally, 
u n ) otherwise. 



[[7.5]] 



where the multipliers u 3 , . . . , i/" that appear in lower order 
5 relaxations used to construct (suboptimal ) solutions of the 
in-dimensional relaxed problem ( [6] 2- [2] 2) have been explicitly 
included. Note that is well-defined since ([6]2.[4]£) can 
always be solved optimally if m=2 . For sufficiently, small 
problems ([6]2. [2]1) or ([6]^. [4]£), one can more efficiently 

10 solve the NP-hard problem by branch and bound. This is the 
reason for the inclusion of the first case; otherwise, the 
relaxed function & 2n is used to guide the non[-] smooth 
optimization phase. That the merit function provides a lower 
bound for the optimal solution follows directly from Fact A. 4 

15 and the following fact. 



Fact A. 5. 



Given the definition of in ([7] 2. [5] 20) 



above j_ 
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^n(u\ • • • / U m /U m+1 7 ...,U n ) [<® m j££nn(u m+1 , • • • / u") f©r all 

multipliers j?, . . . , if, u m+1 , . . . , u n . 
(2.21) 

[[7.6]] 

5 The actual function value used in the optimization phase is 
W^; however, the subgradients are computed as in ([7] 2,. [4119), 
but with the solution i being a suboptimal solution 

constructed from a relaxation procedure applied to the m- 
dimensional problem. Again, the use of these lower order 

10 relaxed problems is the reason for the inclusion of the 
multipliers j?, . . . , if. 

To explain how the merit function is used, suppose there 
is a current multiplier set (uSid/ • • • / u oid) an ^ it is desirable 
to update to a new multiplier set ( ujj+i, . . . , u£ e „) via 

15 (u£i, . . . , uSeJ = Wi, . . . / u n old ) + (Au m+1 , . . . ,Au n ) . Then 
^mn(J£id/ • • • / "Sid/ u old/ - • • / Void) is computed where (u^ id/ . . . , u£ Id ) is 
obtained during the relaxation process used to compute a high 
quality solution to the relaxed /n-dimensional assignment 
problem ([6] 2. [2] 7) at (u m+1 , . . . , u n ) = (u^ d/ - • • / "Sid) • The 

20 decision to accept • • • / uJJ ew ) is then based on 

some other stopping criteria commonly used in line searches. 
Again, (u^ ew , . . . , ujJ e J represents the multiplier set used in the 
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lower level relaxation procedure to construct a high quality- 
feasible solution to the /n-dimensional relaxed problem 
([6]2-[2]7) at (u m+I , . . . f u n ) = (u£i, . . . ,u n new ) . The point is that 
each time one changes {if 1 * 1 , . . . , u n ) and uses the merit function 
5 2mn(i?/ • • • / u™/ u/n+I / - - • / U n ) for comparison purposes, one must 
generally change the lower level multipliers (J?, . . . , i?) . 

An illustration of this merit function for [m=n-l is 
given in "PARTITIONING MULTIPLE DATA SETS, MULTIDIMENSIONAL 
ASSIGNMENTS AND LAGRANGIAN RELAXATION", by Aubrey B. Poore and 

10 Nenad Rijavec, and in "QUADRATIC ASSIGNMENT AND RELATED 
PROBLEMS, DIMACS SERIES IN DISCRETE MATHEMATICS AND 
THEORETICAL COMPUTER SCIENCE", in American Mathematical 
Society, Providence, R.I., Vol. 16, 1994, pp. 25-37 edited by 
Panos 1 m=n-l is given in (A.B. Poore and N. Riiavec. 

15 Partitioning multiple data sets: multidimensional assignments 
and Lagrangian relaxation. In P. M. Pardalos and [Henry 
Wolkowicz . 

] H. Wolkowicz, editors, Oaudratic assignment and related 
problems: DIMACS Series in Discrete Mathematics and 
2 0 Theoretical Computer Science, volume 16, pages 25-37, 1994). 

II. 4. 3 Non- smooth Optimization Methods 
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By Fact A. 4 the function ^(u) is a continuous, piecewise 
affine, and concave, so that the negative of ^(u) is convex. 
Thus the problem of maximizing ^(u) is one of non[-] smooth 
optimization. Since there is a large literature on such 
5 problems, only a brief discussion of the primary classes of 
methods for solving such problems is provided. A first class 
of methods, known as subgradient methods, are reviewed and 
analyzed in [the book by ]JJSLZ. Short, ] ^ Minimization Methods 
for Non-Differentiable Functions[ , ] ^ Springer-Verlag, New 

10 York, 19851. Despite their relative simplicity, subgradient 
methods have some drawbacks that make them inappropriate for 
the tracking problem. They do not guarantee a descent 
direction at each iteration, they lack a clear stopping 
criterion, and exhibit poor convergence (less than linear) . 

15 A more recent and sophisticated class of methods are the 

bundle methods; excellent developments are presented in [the 
books of Hiriart-Urruty and Lemarechal : ]±J.-B. Hiriart-Urruty 
and C. Lemarechal [,] ^ Convex Analysis and Minimization 
Algorithms I and II [ , ] ^ Springer-Verlag, Berlin, 1993; [and 

2 0 the book of ]K.C. Kiwiel [ : 1 . Methods of Descent for Non- 
Differentiable Optimization . In A. Dold and B. Eckmann , 
editors , Lecture Notes in Mathematics, [Vol . 1 volume 1133, 
Berlin, 1985. Springer-Verlag [ , Berlin, 198511- Bundle 
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methods retain a set (or bundle) of previously computed 
subgradients to determine the best possible descent direction 
at the current iteration. Because of the non [-] smoothness of 
the resulting direction may not provide a "sufficient" 
5 decrease in 0^. In this case, bundle algorithms take a "null" 
step, wherein the bundle is enriched by a subgradient close to 
u^. As a result, bundle methods are non-ascent methods which 
satisfy the relaxed descent condition ^(u^) [^mJ <£mn ( u k) i f 
[ u k + i* u J Uk+itMk- These methods have been shown to have good 

10 convergence properties. In particular, bundle method variants 
have been proven to converge in a finite number of steps for 
piecewise affine convex functionals in JJC.C. Kiwiel[, 1 . An 
[A] aggregate [S] subgradient [M] method for [N] non-smooth 
[Cjgonvex [M] minimization [ , 1 . Mathematical Programming^ 27 [, 

15 (1983) pp. 1^320-341, fandl 1983; H. Schramm and J. Zowe[, ] ^ 
A version of the bundle idea for minimizing a non-smooth 
function: Conceptual idea, convergence analysis, numerical 
results!, ] . SI AM Journal on Optimization, 2[ (1992)], [pplNo. 
1:121-152 , February, 1992) . 

20 All of the above non-smooth optimization methods require 

the computation of & m „(u) and a rqe5d> m J a £ d&^ iu) . These in 
turn require the computation of the relaxed cost coefficients 
[ [6]X2. [3] In both the first and second procedures 
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discussed in section II. 3-1, maximization of 0 2n ( u ) must be 
repeatedly evaluated. In the most efficient implementations 
presently known of these two procedures, it was found that at 
least 95% of the computational effort of the entire procedure 
5 is spent in the evaluation of the relaxed cost coefficients 
[ [6] X2. [3] ] JLL as part of computing <P 2n {u) . Thus, generally a 
method should be chosen that makes as efficient use of the 
subgradients and function values as [is ] possible, even at the 
cost sophisticated manipulation of the subgradients. In 

10 evaluating three different bundle procedures: (a) the 
conjugate subgradient method of Wolfe used in section 1.1 of 
the first embodiment of the present invention; (b) the 
aggregate subgradient method of Kiwiel ([described in ]K.C. 
Kiwiel[, ] . A n [A] aggregate [S] j[ubgradient [M] method for 

15 [N] non-smooth [C]gonvex [M] minimization [ , 1 . Mathematical 
Programming^ 27 [, (1983) pp. 1 :320-341 , 1983 ) ; and (c) the 
bundle trust method of Schramm and Zowe ([described in ]H. 
Schramm and J. Zowe[, ] . A version of the bundle idea for 
minimizing a non-smooth function: Conceptual idea, convergence 

20 analysis, numerical results[, ] . SI AM Journal on Optimization, 
2[ (1992)], [pp]gS. 1 : 121-152 , February, 1992 ) , it was 
determined that for a fixed number of non-smooth iterations, 
say, ten, the bundle-trust method provides good quality 
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solutions with the fewest number of function and subgradient 
evaluations of all the methods, and is therefore the preferred 
method. 



5 II. 4. 4 The Two \ 1 ^Dimensional Assignment Problem 

The forward/reverse auction algorithm of Bertsekas and 
Castafion ([as described in ]D.P. Bertsekas and D.A. Castanon[, 
] . A forward/reverse auction algorithm for asymmetric 
assignment problems [, 1 . Computational Optimization and 
10 Applications, 1[ (1992), pp. 1 :277-298 , 1992 ) is used to solve 
the many relaxed two [ ] ^dimensional problems that occur in the 
course of execution. 



II. 4. 5. Initial Multipliers and Hot Starts 
15 The effective use of "hot starts" is fundamental for 

real-time applications. A good initial set of multipliers can 
significantly reduce the number of non-smooth iterations (and 
hence the number of ^ evaluations) required for a high 
quality recovered solution. Further, there are several ways 
20 that multipliers can be reused. First, if an n-dimensional 
problem is relaxed to an zn-dimensional problem, relaxation 
provides the multiplier set {u^ 1 , u m+2 , . . . , u n ) . These can be 
used as the initial multipliers for the [n-m+1 dimensionall (n- 



115 




NUMERI . 0 1 (JSC 2 



m+1) -dimensional recovery problem for n-m+l>2 . This approach 
has also worked well to reduce the number of non-smooth 
iterations during recovery. 

Further, for track maintenance, initial feasible 
5 solutions are generated as follows. When a new scan of 
information (a new observation set) arrives from a sensor, one 
can obtain an initial primal feasible solution by matching new 
reports to existing tracks via a two [ ] ^dimensional assignment 
problem. This is known as the track-while-scan (TWS) 

10 approach. Thus, an initial primal solution exists and then we 
use the above relaxation procedure [is ] to make improvements 
to this TWS solution. Also for track maintenance, multipliers 
are available from a previously solved and closely related 
multi [-] dimensional assignment problems for all but the new 

15 observation set. 



Given a feasible solution of the multi [-] dimensional 
assignment problem, one can consider local search procedures 
20 to improve this result, as described in JC.H. Papadimitriou 
and K. Steiglitz [ , ] ^ Combinatorial Optimization: Algorithms 
and Complex! ty[ ,] ^ Prentice-Hall , Inc . , Englewood Cliffs, NJ, 
1982 [, and] ; J. Pearl [,] ^ Heuristics: [i] intelligent 



II. 4. 6. 



Local Search Methods 
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[s] Search 



[s] Strategies 



for 



[c] Computer 



[p] Problem 



[s] Solvingi ,] ^ Addison-Wesley, Reading, MA, 1984) . One method 
is the idea of Fk- interchanges! k interchanges . Since for 
sparse problems only those active arcs that participate in the 
5 solution are stored, it is difficult to efficiently identify 
eligible variable exchange sets with some data structures for 
solving the assignment problem. However, a local search that 
may be more promising is to use the two [ ] ^dimensional 
assignment solver in the following way. Given a feasible 

10 solution to the multi [-] dimensional assignment problem, the 
indices that correspond to active arcs in the solution are 
enumerated. Subsequently, one coordinate is freed to 
formulate a two- dimensional assignment problem with one index 
corresponding to the enumeration and the other to the freed 

15 coordinate, and then solve a two [ ] ^dimensional assignment 
problem to improve the freed index position. 



2 0 relaxation. Due to the sparsity of assignment problems, 
however, frequently a decomposition of the problem into a 
collection of disjoint components can be done wherein each of 
the components can be solved independently. Due to the setup 



II. 4. 7. 



Problem Decomposition 



The procedures described thus far are all based on 
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costs of Lagrangian relaxation, a branch and bound procedure 
is generally more efficient for small components, say ten to 
twenty feasible 0-1 variables (i.e., z i:L ,,,± n ) . Otherwise, the 
relaxation procedures described above are used. Perhaps the 
5 easiest way to view the decomposition method is to view the 
reports or measurements as a layered graph. A vertex is 
associated with each observation point, and an edge is allowed 
to connect two vertices only if the two observations belong to 
at least one feasible track of observations. Given this 

10 graph, the decomposition problem can then be posed as that of 
identifying the connected subcomponents of a graph which can 
be accomplished by constructing a spanning forest via a depth 
first search algorithm, as described in 1A.V, Aho, J.E. 
Hopcroft, and J.D. Ullmant,]^ The Design and Analysis of 

15 Computer Algorithms[ , ] ^ Addison-Wesley [ Publishing Company, 
Reading] , MA, 19741. Additionally, decomposition of a 

different type might be based on the identification of bi- and 
tri -connected components of a graph and enumerating on the 
connections. Here is a technical explanation. Let 

20 Z = {z ili2 ... fWziIi274 Jz ili2 ... in is not preassigned to zero) 

denote the set of assignable variables. Define an undirected 
graph G{N,A) where the set of nodes is 

N={ [z n |n=l] z?J^=l, . . .,n; i n =l ,N n } 
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and arcs, 

A=( (z^z^) Mn*l] |n*l , j n *0, j^O and there exists z lll2 ... in ez 

such that j n =i n and jf 2 =i 2 } • 
Note that the nodes corresponding to zero index have not 
5 been included in the above defined graph, since two variables 
that have only the zero index in common can be assigned 
independently. Connected components of the graph are then 
easily found by constructing a spanning forest via a depth 
first search. (A detailed algorithm can be found in the book 
10 by Aho, Hopcroft and Ullman cited above) . Furthermore, this 
procedure is used at each level in the relaxation, i.e., is 
applied to each assignment problem [ [3] J^. [1] ] 41 for 
13=3, . . . ,M. 

The original relaxation problem is decomposed first. All 
15 relaxed assignment problems can be decomposed a priori and all 
recovery problems can be decomposed only after they are 
formulated. Hence, in the n-to- (n-1) case, we have n-2 
relaxed problems that can all be decomposed initially, and the 
recovery problems that are not decomposed (since they are all 
20 two [ ] ^dimensional) . In the n-to-2 case, we have only one 
relaxed problem that can be decomposed initially. This case 
yields n-3 recovery problems, which can be decomposed only 
after they are formulated. 
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III. 



NEW APPROACHES TO TRACK INITIATION AND MAINTENANCE 



USING MULTIDIMENSIONAL ASSIGNMENT PROBLEMS 



5 



III.l. 



Introduction 



The ever- increasing demand for sensor surveillance 
systems is to accurate [ly] track and identify objects in real- 
time, even for dense target scenarios and in regions of high 
track contention. The use of multiple sensors, through more 

10 varied information, has the potential to greatly improve 
target state estimation and identification. However, to take 
full advantage of the available data, a multiple frame data 
association approach is needed. The most popular such 
approach is an enumerative technique called multiple 

15 hypothesis tracking (MHT) . As an enumerative technique, MHT 
can be too computationally intensive for real-time needs 
because this problem is NP-hard. A promising approach is to 
utilize the recently developed Lagrangian relaxation 
algorithms (K. P . Pattipati , S. Deb, and Y . Bar-Shalom. A 

2 0 multisensor-multitarget data association algorithm for 
heterogeneous sensors. IEEE Transactions on Aerospace and 
Electronic Systems, 29, No. 2:560-568, April 1993; A.B. Poore 
and N.Riiavec. Partitioning multiple data sets: 
multidimensional assignments and lagrangian relaxation, in 

25 quadratic assignment and related problems, P.M. Pardalos and 
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H.Wolkowicz, editors , DIMACS series in Discrete Mathematics 
and Theoretical Computer Science, 16:25-37, 1994; A.J. 
Robertson III. A class of lagrangian relaxation algorithms for 
the multidimensional problem. Ph.D. Thesis, Colorado State 
5 University, Ft. Collins, CO, 1995 ) for the multidimensional 
assignment problem; however, there are many other potentially 
good approaches to these assignment problems such as LP 
relaxation combined with an interior point method, GRASP, and 
parallel izat ion . 

10 These data association problems are fundamentally 

important in tracking. Thus, our first objective is to 
present an overview of the tracking problem and the context 
within which these data association problems fit. Following 
a brief review of the probabilistic framework for data 

15 association, the second objective is to formulate two models 
for track initiation and maintenance as multidimensional 
assignment problems. This formulation is based on a window 
moving over the frames of data. The first and simpler model 
uses the same window length for track initiation and 

2 0 maintenance, while the second model uses different lengths. 
As the window moves over the frames of data, one creates a 
sequence of these multidimensional assignment problems and 
there is an overlap region of frames common to adjacent 
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windows. Thus, given the solution of one problem, one can 
develop a warm start for the solution to the next problem in 
the sequence, as shown hereinafter. Such information is 
critically important to the design of real-time algorithms. 
5 [ ] 

1 1 1. 2. Overview of the Tracking Problem 
Tracking and data fusion are complex subjects of 
specialization that can only be briefly summarized as they are 
related to the subject of this paper. The processing of track 

10 multiple targets using data from one or more sensors is 
typically partitioned into two stages or major functional 
blocks, namely, signal processing and data [process! processing 
( Y . Bar- Shalom . Mul ti targe t -Mul ti sensor Tracking : Advanced 
Applications. Artech House, Dedham, MA, 1990; Y . Bar- Shalom and 

15 T.E. Fortmann. Tracking and Data Association. Academic Press, 
Boston, MA, 1988; S.S. Blackman. Multiple Target Tracking with 
Radar Applications . Artech House, Dedham, MA, 1986) . The 
first stage is the signal processing that takes raw sensor 
data and outputs reports. Reports are sometimes called 

2 0 observations, threshold exceedances, plots, hits, or returns, 
depending on the type of sensor. The true source of each 
report is usually unknown and can be due to a target of 
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interest, a false signal, or persistent background objects 
that can be moving in the field of view of the sensor. 

Two principal functions of data processing are tracking 
and discrimination or target identification. The 
5 discrimination or target identification function distinguishes 
between tracks that are for targets of interest and other 
tracks such as those due to background. Also, if there is 
enough information, each track is classified as to the type of 
target it appears to represent. Target identification and 

10 discrimination will not be discussed further in this paper 
except to comment here on the use of attributes (including 
identification information) in the data association. 

There are two types of attributes or features, namely, 
discrete valued variables and continuous valued variables. 

15 Available attributes of either or both types should be used in 
computing the log likelihoods or scores for data association. 
In discussing tracking in this paper, the attributes will not 
be mentioned explicitly but it should be understood that some 
reports and some tracks may include attributes information 

20 that should be and can be used in the track processing (both 
data association and filtering) if it is useful. 

[ ] 

III .2.1. Tracking Functions 
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The tracking function accepts reports for each frame of 
data and constructs or maintains a target state estimate, the 
variance -covariance matrix of the state estimation error, and 
refined estimates (or probabilities) of the target attributes. 
5 The state estimate typically includes estimates of target 
position and velocity in three dimensions and possibly also 
acceleration and other variables of interest. 

The tracking function is typically viewed in two stages, 
namely, data association and filtering. Also, the process of 
10 constructing new tracks, called track initiation, is different 
from the process of updating existing tracks, called track 
maintenance . 

In track maintenance, the data association function 
decides how to relate the reports from the current frame of 

15 data to the previously computed tracks. In one approach, at 
most one report is assigned to each track, and in other 
approaches, weights are assigned to the pairings of reports to 
a track. After the data association, the filter function 
updates each target state estimate using the one or more (with 

20 weights) reports that were determined by the data association 
function. A filter commonly used for multiple target tracking 
and data fusion is the well known Kalman filter (or the 
extended Kalman filter in the case of nonlinear problems) or 
a simplified version of it (Y. Bar-Shalom. Multitaraet- 
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Multisensor Tracking: Advanced Applications . Artech House, MA, 
1990; Y. Bar-Shalom and T.E. Fortmann. Tracking and Data 
Association. Academic Press, Boston, MA, 1988; S.S. Blackman. 
Multiple Target Tracking with Radar Applications . Artech 
5 House, Dedham, MA, 1986) . 

In track initiation, typically a sequence of reports is 
selected (one from each of a few frames of data) to construct 
a new track. In track initiation, the filtering function 
constructs a target state estimate and related information 

10 based on the selected sequence of reports. The new track is 
later updated by the track maintenance processing of the 
subsequent frames of reports. 

In some trackers, there is also a tentative tracking 
function for processing recently established tracks until 

15 there is enough confidence to include them in track 
maintenance. While for simplicity of presentation, tentative 
tracking will not be included in the processing discussed in 
this paper, the techniques discussed could readily include one 
or more tentative tracking functions. 

2 0 There have been numerous approaches developed to perform 

the data association function. Since optimal data association 
is far too complex to implement, good but practical sub- 
optimal approaches are pursued. Data association approaches 
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can be classified in a number of ways. One way to classify 
data association approaches is based on the number of data 
frames used in the association process I". 1 (O.E. Drummond. 
Multiple sensor tracking with multiple frame, probabilistic 
5 data association. In Signal and Data Processing of Small 
Targets, SPIE Proceedings, volume 2561, pages 322-336, 1995). 
== In single frame data association for track maintenance, 
"hard decisions" are made [each frame as to how thel on the 
assignment of reports [ are to be related] to the tracks. Some 

10 single frame approaches include: individual nearest neighbor, 
PDA, JPDA, and global nearest neighbor (sequential most 
probable hypothesis tracking) , which uses a two [ ] ^dimensional 
assignment algorith m (Y. Bar-Shalom. Multitarget-Multisensor 
Tracking: Advanced Applications. Artech House, MA, 1990; Y. 

15 Bar-Shalom and T.E. Fortmann. Tracking and Data Association. 
Academic Press, Boston, MA, 1988; S.S. Blackman. Multiple 
Target Tracking with Radar Applications . Artech House, Dedham, 
MA, 1986) . In most single frame data association approaches, 
only one track per obj ect is carried forward to be used in 

2 0 processing the next frame of reports. 

Multiple frame data association is more complex and 
frequently involves purposely [carry! carrying forward more 
than one track per target to be used in processing the next 
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frame of reports. By retaining more than one track per 
target, the tracking performance is improved at the cost of 
increased processing load. The best known multiple frame data 
association approach is an enumerative technique called 
5 multiple hypothesis tracking (MHT) which enumerates all the 
global hypothesis with various pruning rules. 

As shown hereinafter, the log likelihood of the 
probability of a global hypothesis can be decomposed into the 
sum of scores of each track that is contained in the 
10 hypothesis. As a consequence, the most probable hypothesis 
can be identified with the aid of a non- enumerative algorithm 
for the assignment so formulated. 

[ ] III. 2. 2. Interface Between Track Initiation 

And Track Maintenance 

15 There are a number of ways that track initiation and 

track maintenance functions can interact (Q.E . Drummond . 
Multiple target tracking with multiple frame, probabilistic 
data association. In Signal and Data Proceedings, volume 1954, 
pages 394-408, 1993). Two ways that are especially pertinent 

20 to the methods of this paper will be discussed in this 
section . 

One approach is to treat these two functions 
sequentially, by first assigning reports to tracks and then 
using the information that remains to form new tracks. A 
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better approach is to integrate both processes where 
initiating tracks compete for the new frame of reports on an 
equal basis, i.e., at the same time. 

In the integrated approach, the data association for 
5 track initiation and track maintenance processing are combined 
and conducted simultaneously. One assignment array is created 
that includes the track scores for potential tracks for both 
track initiation and track maintenance. The first dimension 
of the assignment problem includes only all the tracks either 

10 created or updated in the processing of the prior frames of 
reports. Each of the remaining dimensions accommodates one 
frame of reports. 

After the assignment algorithm finds a solution, each of 
the previously established tracks are updated with the report 

15 assigned to it from the second dimension of the assignment 
array. The remaining reports in the second dimension that 
were in the assignment solution are firmly established as 
track initiators. The unassigned reports in the second 
dimension of the assignment array are discarded as false 

20 signals. The processing of the next frame of reports repeats 
this process using these updated tracks and the newly 
identified initiators in the first dimension of the cost array 
for processing the new frame of reports. Details of this 
approach are discussed hereinafter. 
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The integrated approach just discussed, uses the same 
number of frames of reports for both track initiation and 
track maintenance. The goal in using the multi -dimensional 
assignment algorithm is to provide improved performance while 
5 minimizing the amount of processing required. Typically, 
track initiation will benefit from more frames of reports than 
will track maintenance. Thus^ a second approach that 
integrates track initiation and track maintenance is discussed 
hereinafter, wherein the number of frames of reports is not 
10 the same for these two functions. This is a novel approach 
that is introduced in this paper. 

[ ] 

III. 2. 3. Multiple Sensor Processing 
There are many advantages and many ways [to]£f combining 
15 data from multiple sensors. There are also many ways of 

categorizing the different algorithmic architectures for 

processing data from multiple sensors. One approach outlines 

four generic algorithmic architectures (O.E. Drummond. 

Multiple sensor tracking with multiple frame, probabilistic 
2 0 data association. In Signal and Data Processing of Small 

Targets, SPIE Proceedings, volume 2561, pages 322-336, 1995) . 

Two of these generic architectures are especially pertinent to 

this paper and are summarized briefly. 
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In the Centralized Fusion algorithmic architecture, 
reports are combined from the various sensors to form global 
tracks. This algorithmic architecture is also called Central 
Level Tracking, Centralized Algorithmic Architecture, or 
5 simply the Type IV algorithmic architecture. In track 
maintenance^ for example, if a single frame data association 
is used, then for data association a frame of reports from one 
sensor is processed with the latest global tracks; then the 
global tracks are updated by the filter function. After the 
10 processing of this frame of reports is completed, a frame of 
the reports from another (or the same) sensor is processed 
with these updated global tracks. This process is continued 
as new frames of data become available to the system as a 
whole . 

15 In Centralized Fusion, using the multi-dimensional 

assignment algorithm for the data association with multiple 
sensors is similar to the processing of data from a single 
sensor. Instead of using multiple frames of reports from a 
single sensor, the multiple frames of reports come from 

20 multiple sensors. The frames of reports from all the sensors 
are ordered based on the nominal time each frame of reports is 
obtained and independent of which sensor provided the reports. 
In this way the approaches discussed in this paper can be 
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extended to process reports from multiple sensors using the 
Centralized Fusion algorithmic architecture. 

The second pertinent algorithmic architecture is Track 
Fusion. This approach is also called the Hierarchical or 
5 Federated algorithmic architecture, sensor level tracking or 
simply the Type II algorithmic architecture. In track 
maintenance^ for example, a processor for each sensor computes 
single sensor tracks; these tracks are then forwarded to the 
global tracker to compute global tracks based on data from all 
10 sensors. 

After the first time a sensor processor forwards tracks 
to the global tracker, then subsequent tracks for the same 
targets are cross-correlated with the existing global tracks. 
This track-to-track cross-correlation is due to the common 

15 history of the current sensor tracks and the tracks from the 
same sensor that were forwarded earlier to the global tracker. 
The processing must take this cross-correlation into account 
and there are a number of ways of compensating for this- cross - 
correlations. One method for dealing with this cross- 

2 0 correlation is to decorrelate the sensor tracks that are sent 
to the global tracker. There are a variety of ways to achieve 
this decorrelation and some are summarized in [Drummond, 1996, 
Signal Data Proc] recent paper (O.E. Drummond. Feedback in 
track fusion without process noise. In Sivnal and Data 
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Processing of Small Targets [ 19951 , SPIE Proceedings , [-2561:] 
volume 2561, pages 369-383 , 1995) . 

Once the sensor tracks are decorrelated they can be 
processed by the global tracker in almost the same way as 



association process is referred to as track (or track-to- 
track) association rather than data (or report-to- track) 



global tracker can process the tracks from the various sensor 
10 processors in much the same way that the global tracker of 
Centralized Fusion processes reports. Accordingly the methods 
described in this paper can be readily extended to processing 
data from multiple sensors using either the Centralized Fusion 
or Track Fusion or even a hybrid combination of both 
15 algorithmic architectures. 

II 1. 3. Formulation of the Data Association Problem 
The goal of this section is to briefly outline the 
probabilistic framework for the data association problems 
2 0 presented in this work. The technical details are presented 
elsewhere (A.B. Poore . Multidimensional assignment formulation 
of data association problems arising from multitarget tracking 
and multisensor data fusion. Computational Optimization and 



5 reports are processed. 



In the case of track fusion the 



association . 



If the sensor tracks are decorrelated, the 
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Applications, 3:27-57, 1994) . The data association problems 
for multisensor and multitarget tracking considered in this 
work are generally posed (A. B . Poore . Multidimensional 
assignment formulation of data association problems arising 
5 from multitarget tracking and multisensor data fusion. 
Computational Optimization and Applications , 3:27-57, 1994) as 
that of maximizing the posterior probability of the 
surveillance region (given the data) according to 



10 where Z N+1 represents N+l data sets, y is a partition of indices 
of the data (and thus induces a partition of the data) T* is 
the finite collection of all such partitions, T is a discrete 
random element defined on r*, [y° is a reference partition, 
] and P(r=Y|2 N+1 ) is the posterior probability of a partition y 

15 being true given the data Z N+1 . The term partition is defined 
below; however, this framework is currently sufficiently 
general to cover set packings and coverings . 

Consider N+l data sets Z(k) (k=l, . . . ,N+1) , each consisting 
of M k actual reports and a dummy report , and let denote the 

2 0 cumulative data set defined by 



Maximize {p(r=y\ Z N+1 ) \y 6 r*}, 



(3.1) 




(3.2) 
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respectively. (The dummy report serves several purposes in 
the representation of missing data, false reports, initiating 
tracks, and terminating tracks (A.B. Poore . Multidimensional 
assignment formulation of data association problems arising 
5 from multitarget tracking and multisensor data fusion. 
Computational Optimization and Applications , 3:27-57, 1994) . ) 
In multisensor data fusion and multitarget tracking the data 
sets Z(k) may represent different classes of objects, and each 
data set can arise from different sensors. For track 

10 initiation the objects are reports that must be partitioned 
into tracks and false alarms. In our formulation of track 
maintenance, which uses a moving window, one data set will be 
tracks and remaining data sets will be reports which are 
assigned to existing tracks, as false reports, or to 

15 initiating tracks. 

We specialize the problem to the case of set partitioning 
(A.B. Poore. Multidimensional assignment formulation of data 
association problems arising from multitarget tracking and 
multisensor data fusion. Computational Optimization and 

2 0 Applications, 3:21-51, 1994) defined in the following way. 
Define a "track of data" as where each i [and ] can assume 
[the] zero or nonzero values [ of 0 and , respectively]. A 
partition of the data will refer to a collection of tracks of 
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data wherein each report occurs exactly once in one of the 
tracks of data and such that all data is used up; the 
occurrence of dummy report is unrestricted. The reference 
partition v° is that in which all reports are declared to be 
5 false. 

Next, under appropriate independence assumptions one can 
show (A.B. Poore . Multidimensional assignment formulation of 
data association problems arising from multitarget tracking 
and multisensor data fusion. Computational Optimization and 
10 Applications, 3 :27-57, 1994) 

P(r = v|z N+1 ) = n L 

p(r = Y °|z N+1 ) tVWev ^--W 

15 L i r .i N+1 ^ s a likelihood ratio containing probabilities for 
detection, maneuvers, and termination as well as probability 
density functions for report errors, track initiation and 
termination . Define 



2 0 c. .... = -InL , 

J [1 if]l f if(z. .... ) are assigned to as a track, 
z, .... = \ i i (3.4) 

x i x w + i l [0 otherwise] 0, otherwise. 



Then, recognizing that the problem (3.1) can be expressed as 
25 the [N+l-dimensional] (N+1) -dimensional assignment problem^ 
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(3.5) ] 



Subject To • ll z ± ± =1 ' i 1 = l/...fM L r (3.5) 

2^ • • • 2^ 2^ • • • 2^ z i 1 „.i 1 

V° i*-i=owo i w+1 =o 1 » +1 
for i k = l,.„,M k and k = 2,.„,N, 

2^ • • - 2^ Z i v ..i„ +1 = 1/ ^W+l = lf 

z . e{0,l} for all i x , 



where c a „ 0 is arbitrarily defined to be zero. Here, each group 
of sums in the constraints represents the fact that each non- 
dummy report occurs exactly once in a "track of data. 1 ' One 
5 can modify this formulation to include mult ^assignments of 
one, some, or all the actual reports. The assignment problem 
(3.5) is changed accordingly. For example, if is to be 
assigned no more than, exactly, or no less than times, then 
the " = 1" in the constraint (3.5) is changed to =, £ , " 

10 respectively. In making these changes, one must pay careful 
attention to the independence assumptions, which need not be 
valid in many applications. Expressions for the likelihood 
ratios can be found in the work of [Poore, 1994, Computational 
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Optimization & Appl . , 3:27-57,1 (A.B. Poore . Multidimensional 
assignments and multitarget tracking, partitioning data sets, 
I.J. Cox, P. Hansen, and B. Jules z, editors. DIMACS Series in 
Discrete Mathematics and Theoretical Computer Science, 
5 American Mathematical Society, Providence, R.I., 19:169-198, 
1995) and the references therein. 

[] 1 1 1. 4. Track Initiation and Maintenance 
In this section we explain two multiframe assignment 

10 formulations to the track initiation and maintenance problem. 
The continued use of all prior information is computationally 
intensive for tracking, so that a window sliding over the 
frames of reports is used as the framework for track 
maintenance and track initiation within the window. The 

15 obj ectives are to describe an improved version of a simple 
method and then to put this into a more general framework in 
which track initiation and maintenance have different length 
moving windows . 
[ ] 

20 

I I I. 4.1. The First Approach to Track 
Maintenance and Initiation 
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The first method as explained in this section is an 
improved version of our first track maintenance scheme (A.B. 
Poore . Multidimensional assignment formulation of data 
association problems arising from multitarget tracking and 
5 multisensor data fusion. Computational Optimization and 
Applications, 3:21-51, 1994) and uses the same window length 
for track initiation and maintenance after the initialization 
step. The process is to start with a window of length N+l 
anchored at frame one. In the first step [one] there is only 

10 one track initiation in that we assume no prior existing 
tracks. In the second and all subsequent frames, there is a 
window of length N anchored at frame k plus a collection of 
tracks up to frame k. This window is denoted by {k[; ]± 
&+1,..., k+N} . The following explanation of the steps is much 

15 like mathematical induction in that we explain the first step 
and then step k to step A:+l . 

[ 1 Track Maintenance and Initiation: Step 1. Let 

{i, ( 1 2 ) , i 2 ( I 2 ) ... , i w ( 1 2 ) , i w+1 ( 1 2 ) ( [ 4 ] 2 . [ 1 ] g) 

be an enumeration of all those zero -one variables in the 
20 solution of the assignment problem (3.5) (i.e., ) excluding 
all the false reports in the solution (i.e., all those zero- 
one variables with exactly one nonzero index) and zero-one 
variables in the solution for which (i lf i 2 ) = (0 , 0) . (The latter 
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can correspond to tracks that initiate on frames three and 
higher.) These denote our initial tracks. 

Consider only the first two index sets in this 
enumeration ( [4] ^. [1] §) and add the zero index I 2 =0 with the 
5 corresponding values of i x and i 2 being zero. Thus, the 
enumeration is now . The notation T 0 (l 0 ) = (z. , 7 z. , 7 ,) will 
be used for this pairing. Suppose now that the next data set, 
i.e., the (N+2) th set, is added to the problem. 

To explain the costs for the new problem, one starts with 
10 the hypothesis that a partition yer* being true is now 
conditioned on the truth of the pairings on the first two 
frames being correct . Corresponding to the sequence 

1 T 0 ( IJ , z . , z . f, the likelihood function is then given by 



15 



20 



L = II L , where 

|r 2 (i 2 ),z i3 ,...,z. w+2 je Y ( [4 ]2. [2a] 7) 

Next, define the cost and the corresponding zero-one variable 
c i i -i = " lni u ...i / 

^2^3 ^+2 ■ L 2 i 3 i W+2 

( [4]3. [2b]£) 

z , ... , ^ | is assigned to Track T 2 (l 2 ) , 



2 3 * 2 ( 0, otherwise, 
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respectively. Then the track maintenance problem Maximize 
{ fLv|ver*1 L r | y£r* ) can be formulated as the following 
multidimensional assignment 



Minimize 



[ ] 



2^3 i H+2 2^3 N+2 



_Subject to: ^ 

2^ ... 2^ Z l,i,...i„ t , = 1 ' - Z 2 = 1 ' ^2' ^ -' L 2' 

1 2 =0 i 4 =0 i Nt2 =0 3 3 »* 2 

2^ 2^ • • • 2^ 2^ • • • .2^. z i,i, ...i„_ 



1 2 =0 i 3 =0 i p . 1= 0 i pn =0 i N+2 =0 2 3 "* 2 

for i p =l,...,M p and p=4, N+2-1, 



E ■■• E z u U ...i= 1 ' i m 2 ^ 1 '-' M - 



W+2 



z n-i 6 {°^ 1 } for all 2 2 ,i 3 ,-,i 



N+2 



(3.9) 



10 Track Maintenance and Initiation: Step k. At the 

beginning of the k th step, we solve the following (N+l) - 
dimensional assignment problem. 

_Minimize[:]_ 22 z2 • • • 22 c i ± ± z i ± -^ww (3 • 1Q) 
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5 



= Subject to:~ 



X > • • • 7 j Z i 4 i = lf V = X ' * * • I*lr? 

x k+l~ u X k+N~ u 

Lv M w , M w „ 

2^ 2^ • • • JL ^i, L H =1 ' -**+l = 1 ' 



x k u x k+2 u X k+N u 



2^ 2^ ••• 2^ i.^ ••• 2-f z i„i I , , . . .1,, „ 

1 =0 i =0 i =0 i =0 i =0 * * +1 * +w 

x k u 1 k+l u x p-l u x p+l u X k+N u 

for i =1,...,/* and p=£+2, W+ir-1, 

{4. [4)] 

2^ 2^ • • • 2^ Z I v i wl -i WH =1 ' ■ i *+N =1 '"-' M k+N' 

L =0 i =0 i =0 * * 1 * N 

L k u x k+l u x k+l+N-2 u 



where for the first step l x and L x are replaced by i x and M lt 
respectively. The first index 2 k in the subscripts corresponds 
to the sequence of tracks where 

W ^^i^ 1 *)^"' 2 !^ 1 *)^ = { z it 1 jt>''"' 2 i^ 1 Jt)} is a track of 
10 data from the solution of the problem prior to the formulation 

of ([4] 3. [4]i0) [ or il . I f k=l db then it is just the first data 

set in frame one. 

Basic Assumption . Suppose problem ( [4] ^. [4] lj)) has been 

solved and let the solution, i.e., those zero-one variables 

15 equal to one, be enumerated by 

{(^(^i)^..^^)--^^^!))}^^! (t 4 U- [5] 11) 
with the following exceptions. 
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a. All zero [ ]^one variables for which (l k , i k+1 ) = (0, 0) are 
excluded. Thus, tracks that initiate on frames after the 
(k+l) th are not included in the list. 

b. All zero-one variables whose subscripts have the form I k =0 
and exactly one nonzero index in the remaining indices 
{ i k+1 , i k + N ) are excluded. These correspond to false 
reports on frames rp=k+l] p - , fk+N] k+N . 

c. All variables [{]z 1 • • [)] for which 

(i^,...,!^) = (0,0,...,0) and i k (l k )=0 finl are excluded from 
( [4]JL. [5] 11) . In other words, the reports on the last 
N+l frames in the are all dummy. Any solution with 
this feature corresponds to a terminated track. 
Given the enumeration ( [4] [5] JJL) , one now fixes the 
assignments on the first two index sets in the list 
( [4] 3 . [5] ljj . The zero index I Jk+1 = 0 is added to the 
enumeration to specify [l k {0) , i Jfe+1 (0))= (0, 0) that is used to 
represent false reports and tracks that initiate on frame k+2 
or later, so that the enumeration ( [4] 3^ . [5] ]JJ is 

n ™ {(^(^)^* + i(^ + i))fc-o- 

[(4.6)] 

Then, for l k+1 = l, L Jt+1 , the such track is denoted by 
and the (N+l) -tuple will denote a track T k+1 (l k+1 ) plus a set 
of reports actual or dummy, that are feasible with the track 
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T k+l (l k+1 ) . [ ]The (N+l) -tuple [ ] will denote a track that 
initiates in the sliding window, i.e., on subsequent frames. 
A false report in the sliding window is one with all but one 
non-zero index i p for some p=k+2, Jc+l+N in the (2\T+1) -tuple 

The corresponding hypothesis about a partition y E r* 
being true is now conditioned on the truth of the L k+1 tracks 
existing at the beginning of the N- frame window. (Thus the 
assignments prior to this sliding window are fixed.) The 
likelihood function is given by 



i *+l**+2'"'**+l+N 



,<0) 



= 1. 



[, L 



k+2 / * 
z i 



k+l+N 



^'([4t3?[7£li2) 



Next, define the cost and the zero-one variable by 



7 = S 

1 1 -i 



1, if{<>",<;:;} is signed to 1^(1^), 

0, otherwise, 



(3.13) 



respectively, so that the track extension problem, which was 
originally formulated as Maximize {l Iyer*}, can be expressed 
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as exactly the same multi-dimensional assignment in (3.4) with 
k replaced by Jc+1 . Thus, we do not repeat it here. 

Ill .4.2. The Second Approach to Track 
5 Maintenance and Initiation 

In the first approach the window lengths were the same 
for both track maintenance and initiation. This can be 
inefficient in that one can usually use a shorter window for 
track maintenance than for track initiation. This section 
10 addresses such a formulation. 

The General Step k. To formulate a problem for track 
initiation and maintenance we consider a moving window 
centered as frame k of length J+J+l denoted by 
[k-I f ... , k, ... , k+J] . In this representation, the window length 
15 for track maintenance is J and that for initiation, X+J+l. The 
objective will be to explain the situation at this center and 
then the move to the same length window at center Jc+l denoted 
by r Tk+l-ll k+1 - I , ... , Jc+1 , ... , [k+l+J] Jc+l+J l , i.e., by moving the 
fframel window to the right one frame at time[ to the right] . 
20 The explanation from the first step follows hereinafter. 
The notation for a track of data is 
End Of Moved Text 
T k U*> ={ Z W ' Z W ' Z w}' ( [4 . 8] 3 . 14) 
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where the index l k is used for an enumeration of those 

reports paired together. We also use the notation T Pl k(l)<) 

to denote the sequence of reports belonging to track T k^^ 

but restricted to frames prior to and including p. Thus, 

W E r *,*< 2 *>' ([4.9]3.15) 
r P ,*< 2 *> = { z W''"' z Vv}' 

r * f *u*> - r P,*u*> u {<^>'"'<w'<u*>} for an y p**- 1 - 



Given this notation for the tracks and partition of the 
data in the frames {k-I, k, k+J} , L T Ptk u k ) will denote the 
accumulated likelihood ratio up to and including frame 
10 p(p<k) for a track that is declared as existing on frame k . 
as a solution of the assignment problem. In this notation, 
the likelihood for T Pf k^ 1 k l and that of the association of 




with track T k^k^ is given by 



15 



( [4 .10] 3 . 16) 
L, 

for any p z k-1, 
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The cost for the assignment of to track 

T k {l k ) and the corresponding zero-one variable are given by 



-ln[L m ,1 
= -ln[L 



r *^Jt' *Jt+l '"^Jt+J 



and 



10 



15 



( [4.11J3.17) 



r I 1, if \ Zi +1 , zf +w \ is assigned to track TA1 J , 

J, Jt i Jc+l"' i Jk+j - - J 



0, otherwise, 



respectively. Likewise, for costs associated with the false 
reports on the frames k-I to k and as associated with 
jz£**, j and the corresponding zero-one variables are 
given by: 



^Jc-l *k x k+l l-k+J 1 k'l" 1 k 1 k + l" 1 j*J 

' 1, if {z ,...,z ,...,z } 

P 1 k-I -t-k^k+l 1 k*J 



( [4.12] 3.18) 



are assigned as a track, 
0, otherwise. 



2 0 For the window of frames in the range {k-I , ... , k) , the easiest 
explanation of the partition of the reports is based on the 
definitions 
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i p \ zf belongs to one of the tracks 1 
1 D k~\ listed on frame k, i.e., to one 

(3.19) 

of the T k (l k ) for some l k =l, lj 
F P ,k = ( {1 M p }\T pik )u{0}, 



so that all of those reports not used in the tracks ^(1^) 
on frame p are put into the set of available reports F pfk for 
the formation of tracks over the entire window. Thus, we 



Minimize 22 22 z2 z2 c ± ± ± ± z l ± ± i 

p^k-I i p 7F pfk r^l z>0 ^-r- ^r-^j W V* + r^w 

Subject To £ £ £ <-r..i*Ww =1 (3 - 20) 

{p=k-I,p*q} i €F pik r=* + l i=0 

for i q <= F p k \{0'} and q=k-I,...,k, 

^ E ^ 52 ^.r-Ww-ib/^ 

for i q =l, M q and q= k+1, k+J, 



r=Jt+l i r =0 * * 1 * J 
2 t =l {r=*+l,r*g) i r =0 * * +1 w 9 * 

for q=k+l, k+J, 
z i± i e {0,1} for all l k zl,i k+1 ,...,i k . 



formulate the assignment problem as 
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Basic Assumptions . Suppose problem ( [4-14] 1 .20 ) has 
been solved and let the solution, i.e. those zero-one 
variables equal to one, be enumerated by 



( [4.1513.21) 

with the following exceptions. 

a. All zero [ ] -one variables in the second list for which 

(i^_j/-f i^f ijt +1 ) = (0^ 0/ 0) are excluded. Thus, tracks 
that initiate on frames after the (Jc+l) th are not 
included in the list. 

b. All false reports a[x]re excluded, i.e., all zero-one 
variables in the second list whose subscripts have 
exactly one nonzero index. 

c. All variables for which = ^' 0) and 
Z{p)nT k (l k ) =0 for p=k-K,...,k where \K^Q] K ^0 is user 
specified. Thus the track T k (l k ) is terminated if it is 
not observed over K+J+l frames. 

Given the enumeration ( [4] 2. [15] 2£) , one now fixes the 
assignments on the all index sets up to and including the 
(k+1) th index sets. 




+i 
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( 



t a ) =< 



T k (l k (l k+1 )),z£(l k+l )) 
for i ^i = 1 ' • • 

for i * + i = 4+i + 1 ' • • -' L * + i- 



(3.22) 



[{4.16}] 

Thus one can now formulate the assignment problem for the 

next problem exactly as in ([4-14J3.20) but with k replaced 

by k+1. Thus, we do not repeat it here. 

5 The Initial Step. Here is one explanation for the 

initial step. First, assume that N=I+J. In this case, we 

start the track initiation with a solution of (3.5) . Let 
{ i x (I 2 ) , i 2 (2 2 ) , i w (I 2 ) , i N+1 (1 2 ) }^ =q 

be an enumeration of the solution set of (3.5), i.e., those 
10 zero-one variables z i 1 a 2 )i 2 (i 2 ):i N+l (i 2 ) =1 , including z 00 ... 0 =l 
corresponding to 1 2 =0, but excluding all those zero-one 
variables that are assigned to one and correspond to false 
reports (i.e., there is exactly one nonzero index in the 
subscript of z i 1 i 2 -i Nn ) , all those zero-one variables that are 
15 assigned to one and correspond to tracks that initiate on 
frames higher than 1+2. Then we fix the data association 
decisions corresponding to the reports in our list of tracks 
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prior to and including frame k+l=I+2. This defines the k for 
problem ([4] 3. [4] 5) and one can then continue the 
development by adding a frame to the window as in the 
general case . 

5 If J+J" > N, then one possibility is to start the 

process with N+l frames, and assuming J<N t proceed as before 

replacing J by N-J for the moment, and continue to add 
frames without lopping off the first frame in the window 
until reaches a window of length T+J+l. Then we proceed as 

10 in the previous paragraph. 

If I+J < N, then one can solve the track initiation 
problem (3.5), formulate the problem with the center of the 
window at k+1 = N+l- J, enumerate the solutions as above, and 
lop off the first N-J-I frames. Then, we proceed just as in 

15 the case I+J=N. 

1 1 1. 5. CONCLUDING COMMENTS 
A primary objective in this work has been to 
demonstrate how multidimensional assignment problems arise 
in the tracking environment. The problem of track 

2 0 initiation and maintenance has been formulated within the 
framework of a moving window over the frames of data. The 
solution of these NP-hard, noisy, large scale, and sparse 
problems to the noise level in the problem is fundamental to 
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superior track estimation and identification. Thus, one 
must utilize the special structure in the problems as well 
take advantage of special information that is available. 
Since these moving windows are overlapping, there are some 
5 algorithm efficiencies that can be identified and that take 
advantages of the overlap in the windows from one frame of 
reports to the next. Here is an example of the use of a 
primal solution of one problem to warm start the solution of 
the next problem in the sequence. 
10 Suppose we have solved problem ([412- [4]]^)) and have 

enumerated all those zero-one variables in the solution of 
([4]2. [4]i£) as in ( [4] 3 . [5] 11) . Add the zero index l k+1 =Q , 
so that the enumeration is 




15 With this enumeration one can define the cost by 




= c 



( [5] 3.24) 



and the two-dimensional assignment problem 



[(5.3)] 



20 Let w be an optimal or feasible solution to this two- 



dimensional assignment problem and define 
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(3.25) 



z 



= < 




and w. . - 1 




5 



( [5]3. [4] 26) 



This need not satisfy the constraints in that there are 
usually many objects left unassigned. Thus, one can 
complete the assignment by using the zero-one variables in 
10 ([4]2. [4]JJ1) with k replaced by k+1 with exactly one nonzero 
index corresponding to any unassigned object or data report. 

For the dual solutions, the multipliers arising from 
the solution of the two-dimensional assignment problem 
( \5 . 31 3 . 25 ) corresponding to the second variable, i.e., 



relaxation scheme (A.B. Poore and N.Riiavec. Partitioning 
multiple data sets: multidimensional assignments and 
Lagrangian relaxation, in quadratic assignment and related 
problems; P.M. Pardalos and H. Wolkowicz, editors, DIMACS 
2 0 series in Discrete Mathematics and Theoretical Computer 
Science, 16:25-37, 1994; A.J. Robertson III. A class of 
lagrangian relaxation algorithms for the multidimensional 
assignment problem. Ph.D. Thesis, Colorado State University, 



15 




These are good initial values to use in a 
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Fr. Collins, CO, 1995) . Finally, note that one can also 
develop a warm start for problem ([4]3. [14J2J)) in a similar 
fashion. 



Multidimensional assignment problems govern the central 
problem of data association in multisensor and multitarget 
tracking, i.e., the problem of partitioning observations 

10 from multiple scans of the surveillance region and from 
either single or multiple sensors into tracks and false 
alarms. This fundamental problem can be stated as (A.B. 
Poore. Multidimensional assignments and multitarget 
tracking, partitioning data sets, I.J. Cox, P. Hansen, and 

15 B. Julesz, editors, DIMACS Series in Discrete Mathematics 
and Theoretical Computer Science, American Mathematical 
Society, Providence, R.I., 19:169-198, 1995; A.B. Poore . 
Multidimensional assignment formulation of data association 
problems arising from multitarget tracking and multisensor 

2 0 data fusion. Computational Optimization and Applications, 
3 :27-57, 1994) 



5 



IV. REVIEW OF NEW RELAXATION SCHEMES 



IV. 1. Introduction 



25 
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Minimize 



1 N ^1 



Subject to: mm (4.1) 

52 . . . 52 z i i =lfi 1 = lf.»/M L f 

i 2 =o i N =o 

• • • 2^ 2^ • • • 2^f z i, i =i 



i 1 = 0 i p _ 1 = 0 i p+1 = 0 i w =0 1 * 

for i p =l,.»,M p and p=2,...,N-l, 



£ ... £ W 1 ' V 1 — '^ 

z. . 6 {0,1} for all i lf ...,i w , 



5 where c 0 ... 0 is arbitrarily defined to be zero and is included 
for notational convenience. One can modify this formulation 
to include mul ti -assignments of one, some, or all of the 
actual reports. The zero index is used in representing 
missing data, false alarms, initiating and terminating 

10 tracks as described in (A.B. Poore . Multidimensional 

assignments and multitarget tracking, partitioning data 
sets, I.J. Cox, P. Hansen, and B. Julesz, editors , DIMACS 
Series in Discrete Mathematics and Theoretical Computer 
Science, American Mathematical Society, Providence, R.I., 

15 19:169-198, 1995; A.B. Poore. Multidimensional assignment 
formulation of data association problems arising from 
multitarget tracking and multisensor data fusion. 
Computational Optimization and Applications, 3:27-57, 1994). 
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In these problems, we assume that all zero-one variables 
z. . with precisely one nonzero index are free to be 

1 1"' 1 N 

assigned and that the corresponding cost coefficients are 
well-defined. (This is a valid assumption in the tracking 
5 environment (A.B. Poore . Multidimensional assignments and 
multitarget tracking, partitioning data sets, I.J. Cox, P. 
Hansen, and B. Julesz, editors, DIMACS Series in Discrete 
Mathematics and Theoretical Computer Science, American 
Mathematical Society, Providence, R.I., 19:169-198, 1995; 

10 A.B. Poore. Multidimensional assignment formulation of data 
association problems arising from multitarget tracking and 
multisensor data fusion. Computational Optimization and 
Applications, 3:27-57, 1994). ) Although not required, these 
cost coefficients with exactly one nonzero index can be 

15 translated to zero by cost shifting (A.B. Poore and N. 
Riiavec. A lagrangian relaxation algorithm for 
multidimensional assignment problems arising form 
multitarget tracking. SI AM Journal of Optimization, 3, No. 
3:544-563, 1993) without changing the optimal assignment. 

20 Finally, this formulation is of sufficient generality to 

include the symmetric problem and the asymmetric inequality 
problem [.] (A.J. Robertson III. A class of lagrangian 
algorithms for the multidimensional assignment problem. 
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Ph.D. Thesis, Colorado State University, Fr. Collins, CO. 
1995) . 

The data association problems in tracking that are 
formulated as Equation ([1]4.1) have several 
characteristics. They are normally sparse, the cost 
coefficients are noisy and the problem is NP-hard (M.R. 
Galey and D.S. Johnson. Computers and Intractability. W.H. 
Freeman and Company, San Francisco, CA, 1979) , but it must 
be solved in real-time. The only known methods for solving 
this NP-hard problem optimally are enumerative in nature, 
with branch- and-bound being the most efficient; however, 
such methods are much too slow for real-time applications. 
Thus one must resort to suboptimal approaches. Ideally, any 
such algorithm should solve the problem to within the noise 
level, assuming, of course, that one can measure this noise 
level in the physical problem and that the mathematical 
method provides a way to decide if the criterion has been 
met . 

There are many algorithms that can be used to construct 
suboptimal solutions to NP-hard combinatorial optimization 
problems. These include greedy (and its many variants), 
relaxation, simulated annealing, tabu search, genetic 



# 



NUMERI .01USC2 

algorithms, and neural network algorithms (C.H. 
Papadimitriou and K.Steiglitz. Combinatorial Optimization: 
Algorithm and Complexity. Prentice-Hall, Inc., Englewood 
Cliffs, NJ f 1982; J. Pearl. Heuristics: Intelligent Search 
5 Strategies for Computer Problem Solving. Addison-Weslev, 
Reading, MA, 1984; C.R. Reeves ed. Modern Heuristic 
Techniques for Combinatorial Problems. Halstead Press, 
Wiley, NY, 1993) . For the three^ dimensional assignment 
problem, Pierskalla (W. Pierskalla. The tri-substitution 

10 method for three-dimensional assignment problem. Journal du 
CORS, 5:71-81, 1967) developed the tri-substitution method, 
which is a variant of the simplex method. Frieze and 
Yadegar (A.M. Frieze and J. Yadegar. An algorithm for 
solving 3 -dimensional assignment problems with application 

15 to scheduling a teaching practice. Journal of the 

Operational Research Society, 39:989-955, 1981) introduced a 
method based on Lagrangian relaxation in which a feasible 
solution is recovered using information provided by the 
relaxed solution. Our choice of approaches is strongly 

2 0 influenced by the need to balance real-time performance and 
solution quality. Lagrangian relaxation based methods have 
been used successfully in prior tracking applications (K.R. 
Pattipati, S. Deb, and Y. Bar-Shalom. A s-dimensional 
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assignment algorithm for track initiation. In Proceedings of 
the IEEE Systems Conference, Kobe, Japan, pages 127-130, 
1992; Y. Bar-Shalom, S. Deb, K.R. Pattipati, and H. 
Tsanakis. A new algorithm for the generalized 
5 multidimensional assignment problem. In Proceedings of the 
IEEE International Conference on Systems, Math, and 
Cybernetics , Chicago, pages 132-136 , 1992; A.B. Poore and N. 
Riiavec. A numerical study of some data association problems 
arising in multitarget tracking. Large Scale Optimization: 

10 State of the Art, W.W. Hager, D.W. Hearn and P.M. Pardalos, 
editors. Kluwer Academic Publishers B.V. , Boston, pages 339- 
361, 1994; A.B. Poore and N. Riiavec. Partitioning multiple 
data sets: multidimensional assignments and lagrangian 
relaxation. In P.M. Pardalos and H. Wolkowicz, editors, 

15 Quadratic assignment and related problems : DIMACS Series in 
Discrete Mathematics and Theoretical Computer Science, . 
volume 16, pages 25-37, 1994; A.B. Poore and N. Riiavec. A 
lagrangian relaxation algorithm for multidimensional 
assignment problems from multitarget tracking. SI AM Journal 

2 0 of Optimization^, No. 3:544-563, 1993). An advantage of 
these methods is that they provide both an upper and lower 
bound on the optimal solution, which can then be used to 
measure the solution quality. These works extend the method 
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of Frieze and Yadegar (A.M. Frieze and J.Yadegar. An 
algorithm for solving 3 -dimensional assignment problems 
with application to scheduling a teaching practice. Journal 
of the Operational Research Society, 32:989-995, 1981) . to 
5 the multidimensional case. 

IV. 2 . Probabilistic Framework for Data Association. 

The goal of this section is to explain the formulation 
10 of the data association problems that governs large classes 
of data association problems in centralized or hybrid 
centralized-sensor level multisensor/multitarget tracking . 
The presentation is brief; technical details are presented 
for both track initiation and maintenance in the work of 
15 [Poore] (A.B. Poore . Multidimensional assignments and 

multitarget tracking, partitioning data sets, I.J. Cox, P. 
Hansen, and B. Julesz, editors, DIMACS Series in Discrete 
Mathematics and Theoretical Computer Science, American 
Mathematical Society, Providence, R.I., 19:169-198, 1995; 
2 0 A.B. Poore. Multidimensional assignment formulation of data 
association problems arising from multitarget tracking and 
multisensor data fusion. Computational Optimization and 
Applications, 3:27-57, 1994) . The formulation is of 
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sufficient generality to cover the MHT work of Reid, 
Blackman and Stein (S.S. Blackman. Multiple Target Tracking 
with Radar Applications . Artech House, Norwood, MA, 1986) 
and modifications by Kurien (T.Kurien. Issues in the 
5 designing of practical multitarget tracking algorithms. In 
Multi target -Multi sensor Tracking: Advanced Applications by 
Y. Bar-Shalom. Artech House, MA, 1990) to include 
maneuvering targets. As suggested by Blackman (S.S. 
Blackman. Multiple Target Tracking with Radar Applications . 

10 Artech House, Norwood, MA, 1986) , this formulation can also 
be modified to include target features (e.g., size and type) 
into the scoring function. The recent work of [Poore and 
Drummond] (A.B. Poore and Q.E. Drummond. Track initiation and 
maintenance using multidimensional assignment problems. In 

15 D.W. Hearn, W.W. Hager, and P.M. Pardalos, editors, Network 
Optimization, volume 450, pages 407-422, Boston, 1996. 
Kluwer Academic Publishers B.V.) significantly extends the 
work of this section to new approaches for multiple sensor 
centralized tracking. Future work will involve extensions 

20 to track-to-track correlation. 

The data association problems for multisensor and 
multitarget tracking considered in this work are generally 
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posed as that of maximizing the posterior probability of the 
surveillance region (given the data) according to 



where Z N represents N data sets, y is a partition of indices 
of the data (and thus induces a partition of the data) , F* 
is the finite collection of all such partitions, T is a 
discrete random element defined on r* # and P{T = y\z N ) is 

10 the posterior probability of a partition y being true given 
the data Z N . The term partition is defined below; however, 
this framework is currently sufficiently general to cover 
set packings and coverings (A.B. Poore . Multidimensional 
assignment formulation of data association problems arising 

15 from multi target tracking and multisensor data fusion. 

Computational Optimization and Applications , 3:21-51, 1994; 

A. B. Poore. Multidimensional assignments and multitarget 
tracking, partitioning data sets, I.J, cox, P. Hansen, and 

B. Julesz, editors, DIMACS Series in Discrete Mathematics 
2 0 and Theoretical Computer Science, American Mathematical 

Society, Providence. R.I., 19:169-198, 1995). 



Maximize 




5 
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Consider N data sets Z(k) (k=l , ... ,N) each of ftk reports 




and let Z N denote the cumulative data set defined 



by 



Z(k) = 




and Z w = {z(l),...,Z(N)}, ( [ 2 ]4. [2] 3) 



5 



respectively. In multisensor data fusion and multitarget 
tracking the data sets Z{k) may represent different classes 
of objects, and each data set can arise from different 
sensors. For track initiation the objects are measurements 

10 that must be partitioned into tracks and false alarms. In 
our formulation of track maintenance (A.B. Poore . 
Multidimensional assignment formulation of data association 
problems arising from multitarget tracking and multisensor 
data fusion. Computational Optimization and Applications , 

15 3:27-57, 1994; A.B. Poore. Multidimensional assignments and 
multitarget tracking, partitioning data sets, I.J, cox, P. 
Hansen, and B, Julesz, editors, DIMACS Series in Discrete 
Mathematics and Theoretical Computer Science, American 
Mathematical Society, Providence, R.I., 19:169-198, 1995) , 

2 0 which uses a moving window over time, one data set will be 
tracks and remaining data sets will be measurements which 
are assigned to existing tracks, as false measurements, or 
to initiating tracks. In sensor level tracking, the objects 
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to be fused are tracks (S.S. Blackman. Multiple Target 
Tracking with Radar Applications . Artech House, Norwood, MA, 
1986) . In centralized fusion (S.S. Blackman. Multiple 
Target Tracking with Radar Applications . Artech House, 
5 Norwood, MA, 1986) , the objects may all be measurements that 
represent targets or false reports, and the problem is to 
determine which measurements emanate from a common source. 



10 partitioning (A.B. Poore . Multidimensional assignment 
formulation of data association problems arising from 
multitarget tracking and multisensor data fusion. 
Computational Optimization and Applications f 3:21-51, 1994; 

A. B. Poore. Multidimensional assignments and multitarget 
15 tracking, partitioning data sets, I.J, cox, P. Hansen, and 

B. Julesz, editors, DIMACS Series in Discrete Mathematics 
and Theoretical Computer Science, American Mathematical 
Society, Providence, R.I., 19:169-198, 1995) defined in the 
following way. First, for notational convenience in 

20 representing tracks, we add a zero index to each of the 

index sets in Equation ([2]4. [2] JO, a dummy report z£ to each 
of the data sets Z(k) in Equation ([2]4. [2]3), and define a 



We specialize the problem to the case of set 



"track of data" as 




where i k can now assume the 
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value [s] of 0 [and respectively], A partition of the data 
will refer to a collection of tracks of data wherein each 
report occurs exactly once in one of the tracks of data and 
such that all data is used up; the occurrence of dummy 
5 report is unrestricted. The dummy report z£ serves several 
purposes in the representation of missing data, false 
reports, initiating tracks, and terminating tracks (A.B. 
Poore . Multidimensional assignment formulation of data 
association problems arising from multitarget tracking and 

10 multisensor data fusion. Computational Optimization and 
Applications, 3:21-51, 1994; A.B. Poore. Multidimensional 
assignments and multitarget tracking, partitioning data 
sets, I.J. cox, P. Hansen, and B. Julesz, editors, DIMACS 
Series in Discrete Mathematics and Theoretical Computer 

15 Science, American Mathematical Society, Providence, R.I,, 

19:169-198, 1995). The reference partition is that in which 
all reports are declared to be false. 

Next under appropriate independence assumptions one can 
2 0 show (A.B. Poore. Multidimensional assignment formulation of 
data association problems arising from multitarget tracking 
and multisensor data fusion. Computational Optimization and 
Applications, 3:21-51, 1994; A.B. Poore. Multidimensional 
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assignments and multitarget tracking, partitioning data 
sets, I.J, cox, P. Hansen, and B. Julesz, editors, DIMACS 
Series in Discrete Mathematics and Theoretical Computer 
Science, American Mathematical Society, Providence, R.I., 
5 19:169-198. 1995) 

P(r= Y °|z N ) Y U^ijey V/ ([2] 4 .[3]4) 

where y° is a reference partition, L i v ..i N is a likelihood 
ratio containing probabilities for detection, maneuvers, and 
10 termination as well as probability density functions for 

measurement errors, track initiation and termination. Then 
-if c. . - - In L . 



15 



-In 



P(r=y|Z N ) 

p ( r = y° I z N ) 



( [2]4. [4] 5) 



Using Equation ([2]£. [3]4)and the zero-one variable 
z i 1 -i N - 1 i-^ (^i/'"^n) 6 Y and 0 otherwise, one can then write 
20 the problem ([2]4.[3]4) as the following N-dimensional 



assignment problem : _ 
[ 



(2.5)] 
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(4.6) 



£ V; 1 (i 2 =l,-,M 2 ), 




p = 2,..., 



w-.i), 



J i i 1"' w 

i 1 J. 2 ...i w _ 1 

z . 6{0,1} for all i lf ...,i 



£ ^,,=1 <vi,...,^>. 



where c 0 ... 0 is arbitrarily defined to be zero. Here, each 
group of sums in the constraints represents the fact that 
each -non -dummy report occurs exactly once in a "track of 
data." One can modify this formulation to include [multi- 
5 assignments] multiassignments of one, some, or all the actual 
reports. The assignment problem Equation ([2]^. [5] 6) is 

k 

changed accordingly. For example, if z± is to be ^assigned 
no more = than, ^exactly, = or = no = less than n ik times, then 
the " = 1" in the constraint ([2]4.[5]<|) is changed to <, =, 
10 ^n^," respectively. Modifications for group tracking and 

multi^resolution features of the surveillance region will be 
addressed in future work. In making these changes, one must 
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pay careful attention to the independence assumptions that 
need not be valid in many applications. 

Expressions for the likelihood ratios i i 1/ ...,i w can be 
found in the work of Poore (A.B. Poore . Multidimensional 
5 assignment formulation of data association problems arising 
from multitarget tracking and multisensor data fusion. 
Computational Optimization and Applications, 3:21-51, 1994; 
A.B. Poore. Multidimensional assignments and multitarget 
tracking, partitioning data sets, I.J, cox, P. Hansen, and 

10 B. Julesz, editors, DIMACS Series in Discrete Mathematics 
and Theoretical Computer Science, American Mathematical 
Society, Providence, R.I., 19:169-198, 1995 ) . These 
expressions include the developments of Reid, Kurien (T. 
Kurien. Issues in the designing of practical multitarget 

15 tracking algorithms. In Multitaraet-Multi sensor Tracking: 
Advanced Applications by Y. Bar-Shalom. Artech House, MA, 
1990) , and Stein and Blackman (S.S. Blackman. Multiple 
Target Tracking with Radar Applications . Artech House, 
Norwood, MA, 1986) . What 1 s more, they are easily modified 

2 0 to include target features and to account for different 
sensor types. In track initiation, the N data sets all 
represent reports from N sensors, possibly all the same. 
For track: maintenance, we use a sliding window of N data 
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sets and one data set containing established tracks. (A.B. 
Poore. Multidimensional assignment formulation of data 
association problems arising from multitarget tracking and 
multisensor data fusion. Computational Optimization and 
5 Applications, 3:27-57, 1994; A.B. Poore. Multidimensional 
assignments and multitarget tracking, partitioning data 
sets, I.J, cox, P. Hansen, and B. Julesz, editors, DIMACS 
Series in Discrete Mathematics and Theoretical Computer 
Science, American Mathematical Society, Providence, R.I., 
10 19 : 169-198 , 1995) . The formulation is the same as above 
except that the dimension of the assignment problem is now 
N+l . [N+l 

To explain the costs for the new problem, one starts 
with the hypothesis that a partition y E F* being true is 
15 now conditioned on the truth of the pairings on the first 
two frames being correct . Corresponding to the sequence , 
the likelihood function is then given by 

Next, define the cost and the corresponding zero-one 
20 variable 

(4.2b) 
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respectively. Then the track maintenance problem Maximize 
{LvlY^r*} can be formulated as the following multidimensional 
] 




10 



£{0,1} for all i^-^i^ 



[ Track Maintence and Initiation: Step k. At the 
beginning of the k th step, we solve] problem 




loop k=l, . . N-2 . At the completion, there remains one 
two-dlmens^ 

2 0 sol^y^ion^^ 



In 
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with the following ( [N+l 5 £jV^&+l) -dimensional assignment 
problem [ . 

Minimize 
Subject to: 

for i p =l, ...,M p and p=k+2 , ... , N, +k+l , (4.4) 



zl k i k+1 -l k+N 6(0,1} for all l k , i k+1 , - , i k+N 

10 

where for the first step l x and L x are replaced by i x and Mj_ 

I-, and L -, are to be replaced by and M , , respectively. [ The 
15 first index l k in the subscripts corresponds to the sequence 

of tracks where is a track of data from the solution of 

the problem prior to the formulation of (4.4) or if k=l then 

it is just the first data set in frame one. 

Suppose problem (4.4) has been solved and let the 
20 solution, i.e., those zero-one variables equal to one, be 

enumerated by 

{ dkdk+i) ,i k+ i(l k+1 ) ,...,i k+N (l k+1 ) ) }L k+1 l k+1 =l (4.5) 

with the following exceptions. 
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a. All zero one variables for which (l k/ i k+1 ) = (0 , 0) are 
excluded. Thus, tracks that initiate on frames after 
the (k+l) th are not included in the list. 

b. All zero-one variables whose subscripts have the form 
5 l k =0 and] 

..... V s N-Jc+l N-*+l 

Minimize J. c i i i z i ± ± 

1 k 1 k+l' ' ' Z N 

Subject To £ i N =l, l k = l,...,L k , (4.7) 

E N-k+l =1 i =1 M 

1 k 1 k+2' * * 1 N 

for i = 1, . . . f Mp and p = Jr+2, . . . , AT-1, 

z w "* +1 =1 i =1 M 



Z I - 



. .6(0,1} for all 1 ,i , . 



roensuret^ 

jyg4££iU*£§yiJ^ exactly one 

nonzero index [in the remaining indices {i k+1/ ... , i^} are 
10 excluded. These correspond to false reports on frames 
p=k+l, ...,k+N. 
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c. All variables z lkik+1 ^ ik+N ) for which (i k + ltmi i k+N ) = (0, 0, . . . , 0) 
and i k (l k )=0 in (4.5). In other words, the reports on 
the last N+l frames in the are all dummy. Any 
solution with this feature corresponds to a terminated 
5 track. 

Given the enumeration (4.5), one now fixes the assignments 
on the first two index sets in the list (4.5). The zero 
index l k¥l = 0 is added to the enumeration to specify that is 
used to represent false reports and tracks that initiate on 
10 frame k+2 or later, so that the enumeration (4.5) is now 

(4.6) 



Then, for , the such track is denoted by , and the 
(N+l) -tuple will denote a track plus a set of reports 

15 actual or dummy, that are feasible with the track . The 
(N+l) -tuple will denote a track that initiates in the 
sliding window, i.e., on subsequent frames. A false report 
in the sliding window is one with all but one non-zero index 
i p for some p=k+2 , ... , k+l+N in the (N+l) -tuple 

2 0 Th el (i.e., variables of the form z i~o*o for 1^=1, . . . . 

and zn'Tn « for jL=1, . . . , m , and p=k+: 



f^eet^beass^ corresponding [hypothesis about 

a partition y G r* being true is now conditioned on the 
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truth of the L k+1 tracks existing at the beginning of the N- 
frame window. (Thus the assignments prior to this sliding 
window are fixed.) The likelihood function is given by 

(4.7a) 

5 Next, define the cost and the zero -one variable by 



0 otherwise, 



(4 . lb) 



respectively, so that the track extension problem, which was 
originally formulated as Maximize {L y /yer* } , can be expressed 
10 as exactly the same multi -dimensional] cost coefficients 




sets. I.J. cox, P. Hansen, and B. Julesz, editors, DIMACS 




lS^lJiJB^^ assignment 
[in (3.4) but with k replaced by k+l. Thus, we do not repeat 
it here. 

2 0 The Second Approach to Track Maintenance and Initiation 

In the first approach the window lengths were the same 
for both track maintenance and initiation. This can be 
inefficient in that one can usually use a shorter window for 
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track maintenance than for track initiation. This section 
addresses such a formulation. 

The General Step To formulate a problem for track 

initiation and maintenance we consider a moving window 
5 centered as frame k of length I+J+l denoted by 

[k-I,...,k,...,k+J] . In this representation, the window length 
for track maintenance is J and that for initiation, I+J+l. 
The objective will be to explain the situation at this 
center and then the move to the same length window at center 
10 k+1 denoted by [k+1 - 1 , ... , k+1 , ... , k+l+J] , i.e., by moving the 

frame one frame at time to the right. The explanation from 
the first step follows hereinafter. 

The notation for a track of data is 

(4.8) 

15 where the index l k is used for an enumeration of those 

reports paired together. We also use the notation T pk (l k ) 
to denote the sequence of reports belonging to track T k (l k ) 
but restricted to frames prior to and including p. Thus, 

(4.9) 

2 0 Given this notation for the tracks and partition of the data 
in the frames {k-I , ... , k, ... , k+j} , LT pk (l k ) will denote the 
accumulated likelihood ratio up to and including frame 
p(p^k) for a track that is declared as existing on frame k 
as a solution of the assignment problem. In this notation, 
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the likelihood for T pk (l k ) and that of the association of 
with track is given by 



(4.10) 

5 respectively . 

The cost for the assignment of to track and the 
corresponding zero one variable are given by 

and (4.11) 

10 

respectively. Likewise, for costs associated with the false 
reports on the frames k-I to k and as associated with and 
the corresponding zero-one variable are given by: 
cFi^-ixi^-ik+J = 
15 zFi^.-ifci^-i^j = (4.12) 

For the window of frames in the range [k-I,...,k], the 

easiest explanation of the partition of the reports is based 
on the definitions 

20 



(4.13) 
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so that all of those reports not used in the tracks on 
frame p are put into the set of available reports F p/k for 
the formation of tracks over the entire window. Thus, we 
formulate the assignment problem as 

Minimize 

Subject to: 



Suppose problem (4.14) has been solved and let the 
solution, i.e., those zero-one variables equal to one, be 
enumerated by 

(4.15) 



with the following exceptions. 

a. All zero one variables in the second list for which are 
excluded. Thus, tracks that initiate on frames after 
the (k+l) th are not included in the list. 

b. All false reports are excluded, i.e., all zero-one 
variables in the second list whose subscripts have 
exactly one nonzero index. 
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c. All variables for which and for p=k-K, k where JC20 
is user specified. Thus the track T k (l k ) is terminated 
if it is not observed over K+J+l frames. 

Given the enumeration (4.15), one now fixes the assignments 
5 on the all index sets up to and including the (k+1) th index 

sets . 



10 Thus one can now formulate the assignment problem for 

the next problem exactly as in (4.14) but with k replaced by 
k + 1 . Thus, we do not repeat it here. 

The Initial Step. Here is one explanation for the 
initial step. First, assume that N=I+J. In this case, we 

15 start the track initiation with a solution of (3.5) . Let 
be an enumeration of the solution set of (3.5), i.e., those 
zero-one variables , including z 00 ... 0 = l corresponding to 1 2 = °/ 
but excluding all those zero-one variables that are assigned 
to one and correspond to false reports (i.e., there is 

20 exactly one nonzero index in the subscript of zi 1 i 2 -i w+1 ) , all 
those zero-one variables that are assigned to one and 
correspond to tracks that initiate on frames higher than 
X+2 . Then we fix the data association decisions 
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corresponding to the reports in our list of tracks prior to 
and including frame k+1 = 1+2. This defines the k for 
problem (4.4) and one can then continue the development by 
adding a frame to the window as in the general case. 
5 If I+J>N, then one possibility is to start the process 

with N+l frames, and assuming J<N, proceed as before 

replacing J by N-J for the moment, and continue to add 
frames without lopping off the first frame in the window 
until reaches a window of length J+J+l. Then we proceed as 

10 in the previous paragraph. 

If I+J<N, then one can solve the track initiation 
problem (3.5), formulate the problem with the center of the 
window at k+1 = N+l- J, enumerate the solutions as above, and 
lop off the first N-J-I frames. Then, we proceed just as in 

15 the case I+J=N. 

A primary objective in this work has been to 
demonstrate how multidimensional assignment problems arise 
in the tracking environment. The problem of track 
initiation and maintenance has been formulated within the 

2 0 framework of a moving window over the frames of data. The 
solution of these NP-hard, noisy, large scale, and sparse 
problems to the noise level in the problem is fundamental to 
superior track estimation and identification. Thus, one 
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must utilize the special structure in the problems as well 
take advantage of special information that is available. 
Since these moving windows are overlapping, there are some 
algorithm efficiencies that can been identified and that 
5 take advantages of the overlap in the windows from one frame 
of reports to the next. Here is an example of the use of a 
primal solution of one problem to warm start the solution of 
the next problem in the sequence. 

Suppose we have solved problem (4.4) and have 
10 enumerated all those zero-one variables in the solution of 
(4.4) as in (4.5). Add the zero index l k+1 =0 , so that the 
enumeration is 

(5.1) 



15 With this enumeration one can define the cost by 

(5.2) 

and the two dimensional assignment problem 
Let w be an optimal or feasible solution to this two- 
dimensional assignment problem and define 
20 (5.4) 



This need not satisfy the constraints in that there are 
usually many objects left unassigned. Thus, one can 
complete the assignment by using the zero-one variables in 
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(4.4) with k replaced by k+1 with exactly one nonzero index 
corresponding to any unassigned object or data report. 

For the dual solutions, the multipliers arising from 
the solution of the two dimensional assignment problem (5.3) 
5 corresponding to the second variable, i.e., . These are 
good initial values for use in a relaxation scheme [11,12]. 
Finally, note that one can also develop a warm start for 
problem (4.14) in a similar fashion. 

Multidimensional assignment problems govern the central 

10 problem of data association in multisensor and multitarget 
tracking, i.e., the problem of partitioning observations 
from multiple scans of the surveillance region and from 
either single or multiple sensors into tracks and false 
alarms. This fundamental problem can be stated as 

15 Minimize 

Subject to: (1.1) 

for i k =l,...,M k and k=2,...,N, 

e{0,l} for all i^-.i^ 

20 

where c 0 ... 0 is arbitrarily defined to be zero and is included 
for notational convenience. One can modify this formulation 
to include multiassignments of one, some, or all of the 
actual reports. The zero index is used in representing 
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missing data, false alarms, initiating and terminating 
tracks. In these problems, we assume that all zero-one 
variables zi 2 *"i w with precisely one nonzero index are free to 
be assigned and that the corresponding cost coefficients are 
5 well-defined. (This is a valid assumption in the tracking 
environment.) Although not required, these cost 
coefficients with exactly one nonzero index can be 
translated to zero by cost shifting without changing the 
optimal assignment . Finally, this formulation is of 

10 sufficient generality to include the symmetric problem and 
the asymmetric inequality problem. 

The data association problems in tracking that are 
formulated as (1.1) have several characteristics. They are 
normally sparse, the cost coefficients are noisy and the 

15 problem is NP-hard, but it must be solved in real-time. The 
only known methods for solving this NP-hard problem 
optimally are enumerative in nature, with branch- and-bound 
being the most efficient; however, such methods are much too 
slow for real-time applications. Thus one must resort to 

20 suboptimal approaches. Ideally, any such algorithm should 
solve the problem to within the noise level, assuming, of 
course, that one can measure this noise level in the 
physical problem and the mathematical method provides a way 
to decide if the criterion has been met. 
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There are many algorithms that can be used to construct 
suboptimal solutions to NP-hard combinatorial optimization 
problems. These include greedy (and its many variants), 
relaxation, simulated annealing, tabu search, genetic 
algorithms, and neural netork algorithms. For the three 
dimensional assignment problem, Pierskalla developed the 
tri-substitution method, which is a variant of the simplex 
method. Frieze and Yadegar introduced a method based on 
Lagrangian relaxation in which a feasible solution is 
recovered using information provided by the relaxed 
solution. Our choice of approaches is strongly influenced 
by the need to balance real-time performance and solution 
quality. Lagrangian relaxation based methods have been used 
successfully in prior tracking applications. An advantage 
of these methods is that they provide both an upper and 
lower bound on the optimal solution, which can then be used 
to measure the solution quality. These works extend the 
method of Frieze and Yadegar to the multidimensional case. 
Probabilistic Framework for Data Association. (ABP) 
The goal of this section is to explain the formulation 
of the data association problems that governs large 
classesf ormulation of data association problems in 
centralized or hybrid centralized-sensor level 
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multisensor/multitarget tracking. The presentation is 
brief; technical details are presented for both track 
initiation and maintenance in the work of Poore. The 
formulation is of sufficient generality to cover the MHT 
5 work of Reid, Blackman and Stein and modifications by Kurien 
to include maneuvering targets. As suggested by Blackman, 
this formulation can also be modified to include target 
features (e.g., size and type) into the scoring function. 
The recent work of Poore and Drummond significantly extends 
10 the work of this section to new approaches for multiple 
sensor centralized tracking. Future work will involve 
extensions to 

track- to- track correlation . 

The data association problems for multisensor 
15 an d] ariging^from multitarget tracking [considered in this 

work are generally posed as that of maximizing the posterior 
probability of the surveillance region (given the data) 
according to 

(2.1) 

2 0 where Z N represents N data sets, y is a partition of indices 
of the data (and thus induces a partition of the data) , T* 
is the finite collection of all such partitions, r is a 
discrete random element defined on T*, y° is a reference 
partition, and P(T=y|z n ) is the posterior probability of a 
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partition y being true given the data Z N . The term 
partition is defined below; however, this framework is 
currently sufficiently general to cover set packings and 
coverings . 

5 Consider N data sets Z(k) (k=I,...,N) each of M k reports , 

and let Z N denote the cumulative data set defined by 

(2.2) 

respectively. Inand multisensor data fusion and multitarget 
tracking the data sets Z(k) may represent different classes 

10 of objects, and each data set can arise from different 

sensors. For track initiation the objects are measurements 
that must be partitioned into tracks and false alarms. In 
our formulation of track maintenance, which uses a moving 
window over time, one data set will be tracks and remaining 

15 data sets will be measurements which are assigned to 

existing tracks, as false measurements, or to initiating 
tracks. In sensor level tracking, the objects to be fused 
are tracks. In centralized fusion, the objects may all be 
measurements that represent targets or false reports, and 

2 0 the problem is to determine which measurements emanate from 
a common source . 

We specialize the problem to the case of set 
partitioning defined in the following way. First, for 
notational convenience in representing tracks, we add a zero 
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index to each of the index sets in (2.2), a dummy report z 
to each of the data sets Z (k) in (2.2), and define a "track 
of data" as where @@i k and z k can now assume the values of 0 
and z, respectively. A partition of the data will refer to a 
5 collection of tracks of data wherein each report occurs 

exactly once in one of the tracks of data and such that all 
data is used up; the occurrence of dummy report is 
unrestricted. The dummy report z serves several purposes 
in the representation of missing data, false reports, 
10 initiating tracks, and terminating tracks. The reference 
partition is that in which all reports are declared to be 
false. 

Next under appropriate independence assumptions one can 

show 

15 (2.3) 
Li lt ..i N is a likelihood ratio containing probabilities for 
detection, maneuvers, and termination as well as probability 
density functions for measurement errors, track initiation 
and termination. Then if ci 1 ...i N =-lnLi 1 ...i N , 

20 (2.4) 
Using (2.3) and the zero-one variable zi 1 ...i N =l if (i 1# ...,i N ) 
6 y an< 3 0 otherwise, one can then write the problem (2.1) as 
the following N-dimensional assignment problem: 

(2.5) 

185 



NUMERI .01USC2 

where c 0 ... 0 is arbitrarily defined to be zero. Here, each 
group of sums in the constraints represents the fact that 
each-non-dummy report occurs exactly once in a "track of 
data." One can modify this formulation to include 
5 multi assignments of one, some, or all the actual reports. 
The assignment problem (2.5) is changed accordingly. For 
example, if z k is to be assigned no more than, exactly, or no 
less than n k times, then the " = 1" in the constraint (2.5) is 
changed to <, ^n k , " respectively. Modifications for 

10 group tracking and multiresolution features of the 

surveillance region will be addressed in future work. In 
making these changes, one must pay careful attention to the 
independence assumptions that need not be valid in many 
applications. 

15 Expressions for the likelihood ratios Li 1# ...,i N# can be 

found in the work of Poore. These expressions include the 
developments of Reid, Kurien, and Stein and Blackman. 
What's more, they are easily modified to include target 
features and to account for different sensor types. In 

20 track initiation, the N data sets all represent reports from 
N sensors, possibly all the same. For track maintenance, we 
use a sliding window of N data sets and one data set 
containing established tracks. The formulation is the same 
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as above except that the dimension of the assignment problem 
is now N+l . 

Overview of the Lagrangian Relaxation Algorithm. (ABP) 

Having discussed the N-dimensional assignment problem (1.1), 
5 we now turn to a description of the Lagrangian relaxation 
algorithm. The algorithm will proceed iteratively for a 
loop k=l, . . ., N-2. At the completion, there remains one 
two-dimensional assignment problem that provides the last 
step which yields an optimal (sometimes) or near-optimal 

10 solution to the original N-dimensional assignment problem. 
In step k of this loop (summarized in Section 4) one starts 
with the following (N-k+1) -dimensional assignment problem 
with one change in notation. If k=l, then the index 
notation 1 1 and L x are to be replaced by i 1 and M x , 

15 respectively. 
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To ensure that a feasible solution of (3.1) always exists 
for a sparse problem, all variables with exactly one nonzero 
5 index (i.e., variables of the form for l k , ...,L k and for 
i p =l, . - . , and p=k+12, . . . , N) assumed free to be assigned with 
the corresponding cost coefficients being well-defined. 
This assumption is valid in the tracking environment. 
Subsection. ] 

10 Computational Optimization and Applications, 3:27-57, 
1994) . 

Section IV. 3.1 presents some of the properties 
associated with the Lagrangian relaxation of ([3]^. [1J7J 
based on relaxing the last (N-k) -sets of constraints to a 

15 two-dimensional one. [ ^l^^^^^^^^^$MS£^M^M: new 
approach to the problem of recovering a high quality 
feasible solution of the original (N-Jc+1) -dimensional 
problem given a feasible solution (optimal or suboptimal) of 
the relaxed two-dimensional problem is described 

20 hereinafter. A summary of the relaxation algorithm is 

given [hereinafter ! in Section IV. 3. 3 , [following which] and 
in^^ect^^^^I^J^^^^ we establish the maximization of the 
Lagrangian dual (an important aspect of the relaxation 
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procedure) to be an unconstrained nonsmooth optimization 
problem and then present a method for computing the 
subgradients . 



IV. 3.1 The Lagranqian Relaxed Assiqnment Problem! 



The {N-k+1) -dimensional problem ( [3^. [lj^) has {N- 
k+1) sets of constraints. A (M p +1 ) -dimensional multiplier 
10 vector (i.e., [u p u^f E [Rm p +i] M^i) associated with the fp- 
th]g^ constraint set will be denoted by 

u p = (u<?, u? , . . . , u M p ) T with Uq = 0 being fixed for each p = k + 

p 

2,...,itf and included for notational convenience only. Now, 
the (N-k+1) [dimensional] -dimensional assignment problem 

15 ( [3] 4 : . [1] 7) is relaxed to a two-dimensional assignment 

problem by incorporating {N-k-1) sets of constraints into 
the objective function via the Lagrangian. Although any {N- 
k-1) constraint sets can be relaxed, we choose the last (N- 
k-1) sets of constraints for convenience. The relaxed 

2 0 problem is 
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One of the major steps in the algorithm is the maximization 
of with respect to the multipliers (u k+2 , . ..,u N ) . It 
turns out that is a concave, continuous, and piecewise 
affine function of the multipliers (u k+2 , ... , u N ) , so that the 
5 maximization of _is a problem of nonsmooth optimization. 
Since many of these algorithms] [???] ] require a function 
value and a subgradient of at any required multiplier 
value (u k+2 , . . u N ) , we address this problem in the next 
subsection. We note, however, that there are other ways to 
10 maximize and the next subsection just addresses one such 
method. 

[ i 

^3^jj^^? r _9 P^^i^s^|^^ Lagrangian Relaxed Assignment 
Problem f . ] 

15 For a function evaluation of , we show that an optimal 

(or suboptimal) solution of this relaxed problem ( [3^. [2] 8.) 
can be constructed from that of a two-dimensional assignment 
problem. Then, the nonsmooth characteristics of are 
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addressed, followed by a method for computing the function 
value and a subgradient. 

Evaluation of . Define for each {l k ,i k+1 ) an index 

(jjk + 2/-/ On) =(jk +2 d k fi^i) J w U*/ijc + i) )^*nd a new cost function 
5 c\ ± by 



10 



Given an index pair {l k , i k+1 ) , ( j k+2 f • • • / Jn) need not be unique, 
resulting in the potential generation of several 
15 subgradients ( [3^. 1 HJ^) . 
Then, 
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* E ft <^A^~t E < (4.10) 

Subject To JJ Vwf 1 ' 1 ^ 1 "'-^' 

S Z lvi k+ i =1 ' i Jc+l =1 ' ' • ,M k+l' 

2 Jl/ if zf'^* 1 ^ ^ = 1 for some , ); ,„ ^ „ , , 

Zj t =\ 1 k 1 k+i i k+2~ i N k+2,-,i N " e{0,l) f or all l k , i k+] 

k k+1 [0, otherwise. 



As an aside, two observations are in order. The first is 
that the search procedure needed for the computation of the 
relaxed cost coefficients in ( [3^. [3 ] 9) is the most 
computationally intensive part of the entire relaxation 
5 algorithm. The second is that a feasible solution z N ' k+1 of a 
sparse problem ( [3] 4 = . [lj^) yields a feasible solution z 2 of 
( [3^. [4^0) via the construction 



2 j 1, if Zi'^i mi = 1 for some (i . 2 , - , i ), 
10 [0, otherwise. 

thus, there are generally solutions other than the one 
nonzero index solution. 

The following Theorem [3^.1 gives a method for 
evaluating and states that an optimal solution of 
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( [3] 4. can be computed from that of ( [3] 4 . [4^10) . 

Furthermore, if the solution of either of these two problems 
is e-optimal, then so is the other. The converse is 

contained in Theorem [3] 4. 2. 

^^^^ 

N-k+1 2 •■Fl- • \ - I • • \ 

B V w i w -r B v w ir ^k+i' - - • V " w* + 2' • • • / J N ) 

and {l k , i M )M0, 0); 

"V^W*^ 0 if {i *+2' (J* + 2' • • • > Jy) 

and U k ,i k+1 )*(0,0); <P N _ k+1 

N-k+1 -i . rr N-k+1 . \^ p ~ 

«ooi M -i. = 1 if c oo^ 2 ..i„ + JL 2 uf p * 0; 

** 2 w p=Jt+2 p 

Theorem [314^1. Let [co 23 ^ be a feasible solution to the 
two^ dimensional assignment problem ( [3^4. [4JJLJD) and define 
fgf' k+1 l^- k+1 by 
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^Wm" 1 ." ^Vw lf (i *+2' • • • / i N ) " ( j k+2 ' ' • • ' 0 N ) 

and (l k ,i k+1 )*(0,0) 
and (l k ,i ktl )*(0,0) 



p=/c+2 



^OOi -i = 0 lf C 00i -i + JLr U i > 



0, 



(4.11) 



Then [of' k+1 1 \^' k+1 is a feasible solution of the Lagrangian 
relaxed problem and 

<t>,r,.A\^-^ l^- k+1 ; [^k£^ f ... t [if ^H-^id^/ 
5 /V +2 Ju*+ 2 , . . . , f^/u^) . 

Tf, in addition, [co-, 1 w 2 is optimal for the two ^dimensional 
problem, then [ of =k+1 ] w N=k+1 is an optimal solution of the 
relaxed problem and 

10 Proof. Let [co N " k+1 lw^± and rco 2 7 wf be as in the hypotheses, 
and let 0 N . ktl ( [co^lw"- 

&£; ft/* a /u** 2 . ru»lu?) and 
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denote the objective function values of ( [314 . [2]_&) and 
( l- 3 i4- ^ ' res P ect i ve ly • Direct verification shows that 
Q N-k + i^M^i sa tisfies the constraints in (34.28) and 

For the remainder of the proof, assume that co 2 wf is optimal 
for ( [3] 4 . [4JJL0) . Let [y"- bl ;x ff -'" J satisfy the constraints in 
([3i4.[2ia) and define 



xli =Ei v ^i N ? +1 i -i for (l.,i k +,)*(0,0) 



x 0 2 0 = l if = (0,0) and 



4> 



-00i w ,-i„ + u i„ 



p=/c+2 



*0 for some (i k+2 ,-, i N ) , 



x 0 2 0 = 0 if U^i^) = (0,0) and 

CoM^-i, + £ "I* >° for all d k+2 ,~, i N ) . 



p=K+2 



Note that [X 2 ^ satisfies^ the constraints in ( [3] 4 . [4] 10J 
and 
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( [X 2 ^/ [^Ju^l, ...,[u»)> N . k+1 (of; rf* 2 , ...,xf) = 0 N _ k+1 (of- 



k+1 



for any feasible solution ry N ~ lt * 1 2 g w ~ Jc * 1 of the constraints 

( [34] . [2] 8J . This implies [(0 s - ktl j w"'** 1 is an optimal solution 

of ( [3] 4 . [2] 8) . Next, 



0, 



N-k+1 



10 ( [u k+21 u^f , . . . , [xf ] \F) 

follows immediately from 

( [ [CO— 1 J^; [Ifju") = „_ 

k+i (co 2 V; u k+23 ^ / fiflxf) 

for an optimal solution [co 2 J^ of ([3]4. [4J10J . 
15 [ 1 _With the exception of one equality being converted to 
an inequality, the following theorem is a converse of 
Theorem [3] 4.1. 

Theorem [3^4.2. Let fof- k+1 l Jbe a feasible solution to 

problem ( [3] 4. [2] 8) and define [co 21 ^ by 
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"li = E w iT\ -i for (l k , i kt ,)* (0,0) , 

1 k+2 " Z N 

Wo 0 = l if {l k ,i k+1 ) = (0,0) and 
$ c oo'Cl-i N + £ < * 0 some (i^ , -, i w ) , 

v * 2 " p=Jt+2 p 

w 0 2 0 = 0 if il k ,i k+1 ) = (0,0) and 

* 2 N p=Jfc+2 p 



10 Then [of lv^_ a feasible solution of the problem ( [31 4. [4 1 10 ) 
and 

k+1 ( [w 2 i^; /V+ 2 7u*+ 2 , [[[l^Ll) . 

If, in addition, 
15 fa/'-** 1 7 w"'** 1 is optimal for ( [3]_4. [2}_8) , then [co± 2 rf is an 
optimal solution of ( [314. [4 1 10 ) . 

k+1 ( [co 2 i^,- r^lu"* 2 , [ifju") and^ 

k+ i ( ru k + 2 7u*+ 2 , . . . , [xfiu») = 

2 0 k+1 ( [U k+2 i^f , . . . # [xfliP) . 

Proof: Let roo N ' k+1 7 /'"^ and [co 2 lwf be as in ^the^ hypotheses, 

and let <P N . k+1 ( \o u - k+1 1 w"-** 1 ; W lv** 2 . . . ., fiflx?) and % K _ 
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k+1 ( [co 2 7 wf_ ; /"l/ +2 7u* +2 , . . . , f i flxF) denote the objective function 
values of ( [314. [218) and ( [314 . K 710 ) , respectively. 
Direct verification shows that [co 2 / w 2 ^ satisfies the 
constraints in ( [314 . [4 710 ) and 

( [co N - k+1 7^-^ ; [v** 2 k£^ f [xfkf) > 5* w _ 

k+1 ( [co 2 i^; [ i/ +2 iu^, . . . , f i/lu?) . 

For the remainder of the proof, assume that [co N ' k+1 7 w"'** 1 is 
optimal for ( [3 ±4. [22_§J and construct [co 2 Ivf as above. Using 
Theorem [311-1, construct co w w_JC+1 jN-k + i j_ by replacing [co N - k+1 J fV / - 
10 in ( [3 74 . 5 711 ) with a)*'""** 1 J™" k+1 l. Then, from that theorem 

<P»-*+i ( ^ N ~ k+1 S N ' k+1 l; [u** 2 !^, lifhl) ^ [ 1= [ lfj N . 
k+1 ( [colV; [ ^ju^f, []_..., [if) < <p N _ k+1 (d*-** 1 ; u* +2 , ... f if). 
Optimality of cf' k+1 7u^ 

Qptim alitv of w""** 1 then implies the last inequality is in 
fact an equality. This proves the theorem. 

The Nonsmooth Optimization Problem. Next, we address 
the nonsmooth properties of the function & N _ 
20 u.i ( [\j*+ 2 ]_u k + 2 , . . . , [ flu* ) as explained in the following 
theorem. 
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[Theorem 3 . 4 /Theorem 4.3. (G.L. Nemhauser and L.A. 




in ( [3M. [217) , let [N Smk+1 ]V^. k+1 (z N ' k+1 ) be the 

~N-k+l 

5 [object 1 objective function value of the (N-k+1) -dimensional 
assignment problem in equation ( [3^. [1 17) corresponding to 
a feasible solution z N ~ k+1 of the constraints in ( [31 4 . [1 17) , 
and let 

Vk+i< uk+2 ' ••• ' u "> ^V H . ktl (z ) *V N _ k+1 (z N - k+1 ) .-50-™ be an 
10 optimal solution of ([34^.17^). Then, <P N . 

*+i ( [u** 2 ! u k+2 * • • . / [ if *! u*^ ) is piecewise affine, concave and 
continuous in { [\f 2 lu k * 2 , . . . , [\fl\f] and 



u=(u* +2 ,...,u»)$ N _ k+1 (u k + 2 , . ,u N ) sl^ +1 (z w -* +1 ) zV N _ k+1 (z N - k+1 ) . 
15 



Most of the algorithms for nonsmooth optimization are based 
on generalized gradients, called smbgrradients, given by the 
following definition for a concave function. 



20 Definition [314. [511. At 50> u = (u 



N - k+1]c+2 (, . . .,u N ) the set is 
called a subdifferential of and is defined for the concave 



function 



a«Wi<«> =<* e RMkt2+1 x...x*^iv k+ i<»> -V k+l («)^ {u)hy 

g T (co-u) for all co e R^ +1 x-xR^ 1 ' N ' k+1 
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gP = a P= u P=O a * W -* +l(U) ={gE 

° ° ° g r (w-u) forall M w * +2+1 x.»xR^ +1 }, 



jvfrere are all permanently fixed. (Recall that these 
5 were used for notational conveience only.) A vector 
g 6 d& N _ k+1 is called a subgradient . 



There is a large literature on such problems, g^g^^^ 
(J^-B^lgiria^ 

1993; K.C. Kiwiel . Methods of descent for nondif f erent iable 




Touzgj^^^ - : 12 1-15.2 L _ .FejDruar^ 

1992; N.Z. Shor. Minimization Methods for Non-Differentiable 
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nondifferentiab^ 

3:145-173, 1975)j_and we have tried a variety of methods 
5 including subgradient methods ^^^^^^^^^^^^^^Q 

New York, 1985) and bundle methods (J^-B^ 
and^C^Lemarechal^^^ 

Algorithms I& XI. Springer-Verlag^ Berlin, 1993^ K.C. 

— ~ ~ — ,_ — — — — — — — — — — — -i~r— — — — — — — — — — — — ~ — *^~„ — ~ — 

10 *y^ij|l^^ 

:_im_i z a: t i on . _ _I n_ _Le_c ture_ Not.es_ _in J^thama^t i^s^ 

Of 



these, we have determined that for a fixed number of 
nonsmooth iterations (e.g., twenty), the bundle trust method 
15 of Schramm and Zowe (H.Schramm and J. Zowe . A_ vers_ion_of__the 



feundle^de^ 

idea, convergence analysis, numerical, results ^_ J? TAM_ Journal. 



on^OjDt^ provides 
excellent quality solutions with the fewest number of 
2 0 function and subgradient evaluations, and is therefore . our 
currently recommended method. 

An Algorithm for Computing <P fc+1 and a Subgradient. 
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Most current software for maximizing the concave function 
requires the value of the function and a subgradient at a 
point . Based on the previous two sections, this can be 
summarized as follows. 
5 1. Starting with problem ( [3] 4 . [1^) , form the relaxed 
problem ( [3] 4 . [2] 8) ; 

2. To solve ([3]4. [2]*}) optimally, defined the two^ 
dimensional assignment problem ( [3] 4 . via the 

transformation ( [3] 4. [3^) ; 

10 3. Solve the two-dimensional problem (3] 4. [4] 10) 

^^^^ 

optimally; 

4. Reconstruct the optimal solution, say ft N ~ k+1 , of 
([3]4. [2^8) via equation ( [3^. [5^9) as in Theorem [3^.1; 

5. Compute the function 




[ ] 6. A subgradient is given by 
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k *i gi P 



dul 



,u N ) 



Jt+1 



-1 



for i p = l, 



, Af p a/id p = /r+2, 



(4. 



Several remarks are in order. First, the optimal solution 
of the minimisation problem ([3]^. [2^8) is required before 
one can remove the minimum sign, replace w z by the 
minimizer and differentiate with respect to to obtain a 
5 subgradient as in ( [3] 4 . 1 [1] 7) . If ft N " k+1 is an approximate 

solution of ([3]4.[2]8), then the subgradient and function 

^^^^ ^^^^ 

values are only approximate with [accruacy] accura^ 
depending on that of 

u k+2 , . . . , u H ) , let { d k d k+1 ) , i k+1 d k+1 ) }^-o^« ir " Jc+1 • Although one 
10 can evaluate the sums in ( [3] 4 . 1 [0] 6) and ( [3] 4 . 1 [l^T) in a 
[straightforward] straight^forward manner, another method is 
based on the following observation. Given the multiplier 
vector (u. x0 , . . . , u„) , [,] let be an enumeration of indices 
w^ k+llk ' ilc+1 } of w 2 (or the first 2 indices of constructed in 
15 equation [314. [5L|) ) such that 
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aa^V^+i) ' ^1(^1)) = (OrO) ^for^ ^regardless of whether 
i p >p/ 0 2 0 = l or not. (The latter can only improve the recovered 
feasible solution.) The evaluation of the bracketed 
quantity in ( [3J^4. 1 [1]J7) for a specific 

k+2,...,N is one minus the number of times the index value i D 
appears in^the ( [y-k+1 ] p-Jc+U th position of the ( [N-k+l ]N- 
k+1) [tuple ] -tuple in the list 

10 Finally, we have presented a method for computing one 

subgradient. If ®#n-*+i (u) is the unique optimal solution of 
( [3J_4 . [2J_8) , then 50. (u)={g} is dif f erentiable at u, g is 

— — v N-k+l 

just the gradient at u, and the subdif f erential ^is a 
singleton. If, corresponding to u,^{ft N_k+1 (1) , . . . , ft N_k+1 (m) } 
15 is not unique, then there are finitely many such solutions, 
say 50 { ^' k+1 (ul) , ... , ^' k+1 (m) } with respective subgradients 
(g(l) , . . . ,g(m)}, [ ] the subdif f erential^ (l k , i k+1 dO (u) is the 
convex hull of {g (1 ) , . . . , g (m) J^JJ^JL^Jgoff^ 
rates^of^jsub^ 

20 programminci . Mathematical Programming. 13:329-347, 1977). 



These nonunique solutions arise in two ways. First, if the 
optimal solution of the two [ 1 ^dimensional assignment 
problem is not unique, then one can genera [1 1 te all optimal 
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solutions of the two [ ] ^dimensional assignment problem 
( [3li. [4 110 ) . Corresponding to each of these solutions and 
to each index pair ( j x , i k+21 f • • • / j N ) i n each solution, compute 
the indices ( jft~ k+lk+2 , - . . , j N ) ^in ( [3]4. [3l|) . If these j's 
5 are not unique, then we can enumerate all the possible 

optimal solutions of ( [3.L4. [2\Q) . Given these solutions, 
one can then compute the corresponding subgradients from 
( [3li.l [112) • 

[ l^In most nonsmooth optimization algorithms, one uses 
10 an [epsilon-subdif f erential] e-subdiff erential . 



Definition [314, [xl 2 . At u = (u k+2 , . . . , u N ) the set 

Od e & N _ k+1 {u) is called an [epsilon subdi f f erential ] e- 
subdiff erential of . , and is defined for the concave 
15 function 

5 V k+ i( u ) =tg e E Hwtl x...x^iv w (") -V k+ i< u > ^ (u) [(u)] 
g T (co-u) + e for all co eR^A. . .x l^ 1 } . W ~* +1 

by 

3 p _ u p _ u p _ 0 3 e <Wi ( ") = { g e x ... x M M " +1 1 * ( w) - * < u , , 

g r (t/-u) + e for ail w el *" 2 * . . . x l M » l >, 
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where are all permanently fixed. (Recall that these were 
used for notational convenicence only.) A vector is called 



IV. 3. 3 The Recovery Procedure. 
The next objective is to explain a recovery procedure, i.e., 
given a feasible (optimal or suboptimal) solution rco 2 / w 2 of 

10 ([314. t 4 110) (° r [co N " k+1 lu£^i of ([314- [2JL&) constructed via 
Theorem 34.1), generate a feasible solution z N ~ k+1 of equation 
( [314- [111) which is close to [co 2 1 in a sense to be 
specified. We first assume that no variables in ( [31 4. [1 ] 7_) 
are pre [s] assigned to zero; this assumption will removed 

15 shortly. The difficulty with the solution rco N ' k+1 1 w"- k+1 in 
Theorem 314-1 is that it need not satisfy the last ( [N-k- 
l ]N-k-l ) sets of constraints in ( [314. [117) . (Note, 
however, that if [tf lw 2 is an optimal solution ( [314- f4 7 10 ) 
and [co N " k+1 1 w**~ k+1 , as constructed in Theorem 34^1, satisfies 

20 the relaxed constraints, then [co N ~ k+1 7 w""** 1 is optimal for 
( [314- [11^) •) T ^ e recovery proceudregrocedu^ described 





[HOW DO YOU COMPUTE 



5 




here is designed to preserve the 0-1 character of the 



solution fco 2 ]w^_ of ( [3 /4 . [4 /10 ) as far as possible: If 
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z 2 lk =1 and l k *0 or i N _ k+lk+1 =1 * 0 , the corresponding 

"lo-ipO-o^-k+i 

feasible solution z N ~ k+1 of ( [314. [117) is construed so that 

z oot + k t 21 ---Vi N = 1 ^ for some ^^2/.--/iJ. By this reasoning, 
variables of the form cd 2 */-* + i ^can be assigned a value of 

O00Oijt +2 = l~i N """ 

5 one in the recovery problem only if . However, variables 
co 2 N-jt + i ^will be treated differently in the recovery 

ooooi| +2 ...i w 

procedure in that they can be assigned either zero or one 

L k 1 2 

independent of the value [of] { (l k (l w ) f i k+1 hj^ooo • 

This increases the feasible set of the recovery ^problem, 
10 leading^ to a potentially better solution. 

[ ]Let {U k (l k ^)i i k+1 il k+l )}^\ J>^ an enumeration ^of^ 
indices ^^{Ik* ik + i) of [co 2 (or the first 2 indices of [oo N ~ 
K+1 lw N ' k+1 constructed in equation ( [5) ) such that for and for 
l k+1 =0 74 . 11) ) such that 
15 wf (1 , , {1 ,=1 for^ 



(1.(1.^ ) , i. - (1. . ) ) * (0, 0) and 



Ci i i = Ci / i \ i n h i ror k = i 

2 3 * * " N i * 2' 2 ' 2' 3* " N 

regardless of whether 

N-k N-k+l 

c • ~~ c 

^"k+lAk+2 " " * ^"N 1 k 1 k+l^k+l 1 k+l^-k+2" * " *N 2 - 

w 00 = l or 

N 

20 not. To defihe 2 -t^ i2 j>^^ i3 &^ 

for k * 2 and l k = 0, . . . ,L k 
that resores feasibility, let 
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c V3...i M - Ci i u 2 )i 2 a 2 )i 3 ...i N for) 
c Nk = c Nk+1 

= c N 

i l< 1 2...k+0 i 2t 1 2...k+l) i 3< 1 3.. 



for kz2 and I„=0, . 
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where 

...* + i = V-*-i 0 ■■■ ° V ) = J J ^ ( - ( h < ) ) - > ) (4.19: 



for jn = 2, . . . Jc+1 and where [o] "°" denotes function 
composition. For example^ 1 23 = 2 2 d3) and 2 234 [= l 2 o 1 3 (1 4 )] = 

Then the {N-k) -dimensional assignment problem that 
restores feasibility is_ 

[ (3.14)] 



. . N-k N-k 

Minimize 7. c i i i z i i i 

Subject To £ . Z W M .-^ =1 ' I jfi = 1 '--" I W 

1 Jt+2 • * • 2 W 

E . <;^ 2 ...i N =i. i M =i,... f » 4t2 , (4.20) 

^♦l 2 *-^" " ' 2 N 

EN-k - 
Z l i i =1 

2 /c 2 Jt+l * * " 2 p-l 2 p+l " • * 2 N 

for i = 1, . . W M and p = £+3, . . . , N-l, 

I Jt+l 1 Jt+2* * " 2 W-1 
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; Let Y be an optimal or feasible solution to this [N - 
k .1 (N-k) - dimensional assignment problem. The recovered 
feasible solution z N is defined by 



(3.14)1 



10 



Said in a different way, the recovered f reasible 7 feasible 
solution z N z N corresponding to the multiplier set 



{ \u k + 2 lu k+2 . f... 1 u N ) is then defined by 

N-k+1 



1, if r a _ =1, 

0, otherwise, 



where l m ... k+1 is defined in ([3il.l[3l|) and r° / "°" denotes 
function composition. 

The next objective is to remove the assumption that all 

15 cost coefficients are defined and all zero-one variables are 
free to be assigned. We first note that the above 
construction of a reduced dimension assignment problem 
( [324- I"14 7 20 ) is valid as long as all cost coefficients c* 7 "* 
are defined and all zero-one variables in z N ' k are free to be 

20 assigned. Modifications are necessary for sparse problems. 
If the cost coefficient is well-defined and the zero-one 
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variable is free to be assigned to zero or one, then 
define as in ([314.1 [228) with 

being free to be assigned. Otherwise, is preassigned to 
zero or the corresponding arc is not allowed in the feasible 
5 set of arcs. To ensure that a feasible solution exists, we 
now only need ensure that the variables are free to be 
assigned for ri kfl =0 ll^i =0 .1,..., Lu. , with the corresponding 
cost coefficients being well-defined. (Recall that all 
variables of the form and are free (to be assigned) ^and 

10 all coefficients of the form and [ 2 

AT- Jt+1 N ~ k+l 
Ci / -j \ 7 - /] \n n are well-defined for n h =0 / l k 

^0 i ,...,L k and ri p =0 7i r =0,...,M r for p=k+l , ... ,N.) This is 

accomplished as follows: If the cost coefficient is well- 

N-k+1 N ~ k 
defined and is free, then define c 7 / 7 \ (1 \ n . n 

15 with being free. Otherwise, since all variables of the 

form [ J.and are free to be assigned with the 

corresponding costs being well-defined, set 

C Q>-.0 = ?Q~oQ'"~&^-^ Wl 2 ere ,4f St t8rm iS omitted if 



20 



lkdk + i)=0 an^-Che^cond, if (l k+1 ) =0 . For l k (l k+1 )=0 and 



j^tfo 0 ' d £ il ^ i =1 1 =1 L 



&&&&&&&& 



z]7 ^K 1 } for 311 h-i' V c o""o* = CoT 1 - 
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The description of the algorithm ends with [k=N+2 7Je = 
N+2. L The resulting recovery problem defined in 
([314.1 [41QJ with [k=N-2 Ik = N-2 is 

L l = 1 l ( "^2 ■ N-l ^ ' 1 Z = 1 2 ( " i 2-N-l ^ ' 2 3 = ■ L 3 ^ -^3-M-l ^ ' " ' 

srwise. 

^i^l, 1^ = 1,...^!' ( 4 - 22 > 

_Let Y be an optimal or feasible solution to this [2- 
5 dimensional 7 two-dimensional assignment problem. The 
recovered feasible solution z N is defined by 
[ I 

[N 1, if = l x (-i 2 .-N-l^ ' 2 2 = 1 2 ^2 -W-l^ ' 2 3 = 2 3 (■ i 3-N-l) ' " 

^ -» i W-2 = i N-2( 1 W-2W-l^ i W-l = i W-l^W-l) and ^-W^ 1 

or if (I w . lf i w ) = (0,0) , 



( r3.1874.23) 



0, otherwise. 



Said in a different way, the recovered feasible solution z N 
is then defined by 
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1, if Y x i =1 
0, otherwise. 



[N 

(3.19)1 

where l m ... N . 2 is defined in ([314-1 [3^9) and \° 1 "<*" denotes 
function composition. To complete the description, let be 
5 an enumeration of indices i r |I w _ 1 (I w ) f i N (l w ) i w } = 1 for [ lof 
such that (y i n , . n ,L /n ft . anH =l_for 

\ l N _ 1 ( 1 N ) , i N \l N )) * \0r 0) and 

(l WJ (2J ,i M (l M )) = "(0.0) for I N =0 and 

Vl = ^ Vl N ) ' 1 N= 1* ( 2 N ) '' I Vl ( V ' ^ ( V ) = ( 0 ' 

) otherwise. 

regardless of whether or not . Then the recovered 

10 feasible solution can be written as 



tN} N * N 9 » 

W W Wj™ 1 . 



1 , if i 1 = i, ( I 2mW ) , i 2 = i 2 ( I 2mW ) , i 3 = 
0, otherwise. 



15 IV. 3. 5. The Upper and Lower [Bound. 7 Bounds 

The upper bound on the feasible solution is given by 



( [314.2 [116) 
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and the lower by [$^2^ (u 3 , u N ) for any multiplier value 
(u 3 ,..., fu N l u N ) . In particular, we have 



<P N (u\...,u N ) z <P N (u 3 ,...,u N )l V N (z N )l V N (^ w )0[3i4.2[2i7) 

5 

where {u 3 ,...,u N ) is any multiplier value, (u 3 ,...,^) is a 
maximizer of rd^/jg^Cu 3 , u N ) , 2? is an optimal solution of 
the multidimensional assignment problem ( [324 . ^12) anc ^ ^ 
is the recovered feasible solution defined by ( [3 14 . 2 TO 7_5) . 
10 Finally, we conclude with the observation that V N {z)=V 2 {Y) 
where Y is the optimal solution of ( \3 . 17 74 . 23 ) used in the 
construction of z in equations ( [314 . [18121) - ( [314 .2 [Olg) . 

IV. 3. 6. Reuse of Multipliers f . 7 
15 _ Since the most computationally expensive part of the 
algorithm occurs in the maximization of ^-ic+i' fc he 
development of algorithms that can make use of hot starts or 
multipliers close to the optimal are fundamentally important 
for real-time speed. The purpose of this section is to 
20 demonstrate that the multiplier set obtained at stage kzl 

provide good starting values for those obtained at step Jc+1. 

Theorem [3 74 . [xx74 , Let (u 3 ,...,u N ) denote 
a multiplier set obtained in the f maximu m 7 maximi zation of 

214 



NUMERI .01 USC 2 

0 N [ CJ_(u 3 , u N ) . Then this multiplier set satisfies the 

string of inequalities 

[ ,L 

...,u N ' w -* +1 )0 w (u 3 , . . .,u N ) ^ N .j(u 4 , . . .,u N ) z . . . *<P 3 (u w ) s(P 2 sV ff (z) , 

where after the first step in the maximization of f& N l the 
5 multipliers are not changed in the remaining steps. 
Furthermore, to improve this inequiality, let 
® l u 2 * k ' N ~ k * 1 (u k+2 ,..., u*** 1 '" k l ) . denote a maximizer of 
<V* +1 (u k+2 r-f u N ) • Then, we have 



10 [ (3.xxx) 

Inequalities depend on solution of 2d assignment 
problem. J. 

*»-x*i • • ■ ' " N > * (u 2+ *' N "* +1 , . . . , u N ' N '* +1 ) (4.29) 

£<P N _^(u 3+ *' N -* +1 , ...,u w ' N "* +1 ) 
i* w _ Jk (u 3+ *' M -* / ... f u 
s . . . <; d> 3 (u N ' 3 ) s <i> 2 <: V N (2) . 

Proof: Suppose we have a value of the multipliers 
15 {u k+2 , u N ) obtained in the maximization of 
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*»- k +il uk+2 ' Minimize <p N _ k+1 (z w "* +1 ; u* +2 , . . . , u N ) 

Subject To 2 z?-£\ iM „ mi =l, l k = l,...L k , 

1 k+l 1 k+2* ' ' X N 

EW-Jt+1 . _ - M 



(4.30) 



where 



[(3.xx)2 



<i> (z H -* +1 ;u* +2 f . 

*H-k*l 



EW-Jt + 1 N-Jt+1 



p=k+2 i p =0 



,u w )u w ) 



N-Jfc + 1 

z l i i ~ 1 



N-k+1 , ^ p W-Jt+1 ^ \^ p 



p=Jt+2 



(4.31) 



[ J_These need not be the maximizer; however, we do assume 
that we have solved the above minimization problem optimally 

to evaluate (u t U w ),i w {i w )^ 1 s0 ^ |oW '-» B ' 1 ' JuSt as 
in the definition of the [earlier j_recovery problem 

10 discussed earlier, let be an 
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enumeration of indices { l k , i k+1 [ ) ]_X of the optimal solution w 2 
of problem ( [314.. f4 710 ) (or the first 2 indices of the 
solution \v^ =k+1 7 w"'** 1 of the relaxed problem ( [314 . [2JQ) 
constructed in equation ( [314. t 5 lJ=U ) such that 

(l k U k . 1 ).i k . 1 U k . 1 ))*'H0,0)axia 0 s l(l 4 (l w ),i w (l w )) s IO, for 1^=0 

<t> ( U* +2 , U ) s 

re^atdless of 7 whether 



Minimize ^^{z^ 1 ; u k + 2 , u N ) Qr 

10 Then, the final value of the objective function can be 

expressed as 



where 



N-k+1 A p 

C I (1 )i (1 )i -i + 2^ u i 



15 



//> / _N-*+l. „*+2 , ..Hi . 



7V-A-+1 

Zl *<- Z Jk + l) 1 **l( 1 *-l) i * + 2- i W 

[(3.xxx)l 



If now another constrain set, say the (^+2) th set, is added 
to this problem, one has 
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Since the constraint 

or i^. +2 = l , ...,M^ +2 is now imposed, the feasible region is 
5 smaller, so that one has 

(u* +2 ,...,u N )s <P N _ k+1 {u k+2 ,...,u N ) 

s ^ +1 (u^ 3 ,...,u w ), 

where the fact that 

3> w (u 3 ,...,u w ) £<£ N _i( u4 '-"' uN > *«.** 3 (u w ) ^d> 2 H y w (£ ) # N=Jt+1 ^does not 
10 depend on the multiplier set u k+2 due the added constraint 
set. Also, the last equivalence follows from the fact that 

^®N-k\il ®n-k\i ( uk+3 fu N lv^) is the relaxed problem[ (3.xx)j_ 

(with k replaced by Jc+1) for the recovery problem [ (3.xx)j_ 
in step k of the above algorithm. Continuing this way, one 
15 can compute (u 3 ,...,^) at the first step (k=l) , fix them 

thereafter, and perform no optimization at the subsequent 
steps to obtain 



(u 2+JC ' w -* +1 ,..., ^u^..,u w )^ N . i (u^...,u N )^■•^ 3 (u w )^ 2 ^ w (^) / 
where <P 2 is defined in ( [314. ri7 722 ) . 
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To explain how to improve this inequality, let 

3> ( u 2+/c,w ~* +1 (u k+2 , u nW ' w . denote a maximizer of 

^ w -/c + i ( u k+2 , u w ) . Then by the same reasoning one has 

l"(3.xx) /the result as stated in the theorem. [Q.E.D.J, 

Summary of the Lagrangian Relaxation Problem. (APB)2 

V. SUMMARY OF THE LAGRANGIAN RELAXATION PROBLEM 

Given the multidimensional assignment problem 

[ (4.1)1 



Minimize E ^...^...i. 



1 • * * N 



Subject To E z ,- i =1 (i, = 1, • • -WJ , (5.1) 
i 2 i,...i„ 

E z i,...i„ =1 (i 2 = l,...Mj)/ 



z =1 (i =l,...tf and p = 2, . . . , N-l) , 



J] z. . = 1, (i w =lf . . ) , 

l x l 2 . . . l^J 

z. 6 {0,1} for all i , . . . i , 

i 1# . . . i w « 

10 where is arbitrarily defined to be zero and is included 

for notational convenience and where the superscript N on 

both c and z is not an exponent, but denotes the dimension 

of the subscripts and the assignment problem as stated in 
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([2] 3^5). In relaxed and recovery problems ^need not be 
zero! In this problem, we assume that all zero-one 

z k * 2N 

variables (u ir " iw , u w ) , ^with precisely one nonzero index 
are free to be assigned and that the corresponding cost 
5 coefficients are well-defined. (This is a valid assumption 
in the tracking environment (A.B. Poore . Multidimensional 

sets, I.J. Cox, P. Hansen, and B. Julesz, edit^rSj_DIMACS 



10 




19:169-198, 1995; A.B. Poore. Multidimensional assignment 




Tomgutatio^ _Apnlications J _ _3_:_2 7_-_57^ 

15 ijy^y^* ) Although not required, these cost coefficients with 
exactly one nonzero index can be translated to zero by cost 
shifting (A.B. Poore and N. Riiavec. A lagrangian 



re. 



;ion algorithm tor multidimensional assignment 



problems arising from multitarget trackinc 



Tournal ol 



without changing the 



2 0 2E£^S4^^^SA«tfife«£S:- 
optimal assignment . 

Having explained many of the relaxation features, it is 

now appropriate to present the Lagrangian relaxation 
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algorithm, which iteratively constructs a feasible solution 
to the Itf-dimensional assignment problem ([4] 5..1). 

Algorithm [4] 5.1 (Lagrangian Relaxation Algorithm) [ .] 
To construct a high quality feasible solution, [enoted] 
5 denoted by z N , of the assignment problem ( [4] ^.1) , proceed 
as follows: 

0. Initialize the multipliers^(u k+2 ,...,u N ) = (0,...,0). e.g., 

(u k+2 ,...,u N ) = (0,...,0) . 
For k-l,...,N-2, do 
10 1. For the Lagrangian relaxed problem ([3] ^.[2] 8) from 

the problem ( [3] 4. [1] 7) by relaxing the last (N-k-1) 

^^^^ 

sets of constraints . 
2. Use a nonsmooth optimization technique to solve 

15 [ (4.2)] 



(u* +2 ,...,u w ) Maximize { <p N _ k+1 \ u p et" p+1 for p = k+2,...,N 

with u 0 p =0 being fixed), 

where $ N _ k+1 ( u k+2 , u w ) , is defined by equation ([3]^. [2] 
8). The algorithm [following problem (] in Section 
' provides one method for computing a 

2 0 function value and a subgradient out of the 
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Jt+2 N W 

subdifferential at (Q ), as required in most 

nonsmooth optimization techniques. 

3. Given an approximate or optimal maximizer of ([4] ^.2), 
say (Q k+2 ,-, Q N ), let [w 2 ] <J> 1: L . 1 (u^" u **) denote the 

5 optimal solution of the two-dimensional assignment 

problem ([3] 4. [4] 10) corresponding to this maximizer 

of •* tal V*i<« 3 «* +2 '"'«"> * s V N (Z N ) . 

4. Formulate the recovery {N-k) -dimensional problem ([3] 

4.1 [3] ^9) , modified as discussed at the end of 
10 subsection [ (] ^TV^3 . 3 [)] for sparse problems. At this 

stage, z N , as defined in ([3] 4. [15] 21), contains the 
alignment of the indices (i lf ... , i^ +1 } . 

Formulate the final two-dimensional problem ([3] 4. [17] 22) 

^^^^^ 

and complete the final recovered solution z" as in ( [3] 
15 4.2[0] 5). 

To explain the lower and upper bounds, let <P N be as 
defined in ([3] 4. [2] 8) with k=l t let V N (z N ) be the 
objective function value of the N-dimensional assignment 
problem in equation ([2] 3^5) corresponding to a feasible 
20 solution z N of the constraints in ([2] 3.5), and let ^~ k+1 be 
an optimal solution of ([2] 3.5). Then, 
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?0 $0 C W/W, 

t to 22 }2 z wr 1 ' i i =1 " 

i,=0 i,=0 1 2 3 



' Z iii = 1 ' i 2 = 1 '- ' M 2' 

e (0, 1) for all i 1 i 2 i 3 . 



► z 
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(45.3) 



is the desired inequality. 
COMMENTS : 

5 1. Step 2 is the computational part of the algorithm. 

Evaluating & N _k+i and computing a subgradient used in 
the search procedure requires 99% of the computing time 
in the algorithm. This part uses two ^dimensional 
assignment algorithms, a search over a large number of 
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indices, and a nonsmooth optimization algorithm. It is 
the second part (the search) that consumes 99% of the 
computational time and this is almost entirely 
parallelizable . 

5 2. In track maintenance, the warm starts for the initial 
multipliers for the first step are available. For the 
relaxation procedure, initial multipliers are available 
in step 2 from the prior step in the algorithm. 
3. There are many variations on the above algorithm. For 

10 example, one can compute a good solution on the first 

pass {k=l) and not perform the optimization in (2) 
thereafter. This yields a great solution. Thus one 
can continue the optimization at the first pass, and 
immediately compute quality feasible solutions to the 

15 problem. 

[ Lagrangian Relaxation Algorithm for the 3 -Dimensional 
Algorithm. (ABP) ] VI. LAGRANGIAN RELAXATION 

ALGORITHM.^ 

FOR THE THREE-DIMENSIONAL ALGORITHM 



2 0 In this section, we illustrate the relaxation algorithm 

for the three ^.dimensional assignment problem. Having, 
discussed the general relaxation scheme, 
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, z , , and z 0Q , Minimize £ £ 2 c m z - 
Subject To £j z i 1 i 2 i 3 =1 ' i^l,-,^, 

^ ^ Z W. =1 ' i 2 = 1 '-' W 2' 

1^=0 i 3 =0 1 2 3 

^V,^ 0 ' 1 ) for all i,^!,. 



( [5]6.1) 



To ensure that a feasible solution of ([5]£.l) always exists 
5 for a sparse problem, all variables with exactly one nonzero 
index (i.e., variables of the form for i D =l,...,M D and 

^^^^ ir ir 

p-1,2,3 are assumed free to be assigned with the 
corresponding cost coefficients being well -defined^. This 
assumption is valid in the tracking environment [ [**refPoa, 
10 **ref Pob] . ] 

[ 1 1 A ._B_._ _ Poore ^ _ Mul t idimens_ional. _a_ss iqnment s_ _anc 




Ma£hemat^ic^^ 



15 
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5 The three -dimensional assignment problem ([5] 6.1) has 

three sets of constraints and one can describe the relaxation 
by relaxing any of the three sets of constraints, the 
description here is based on relaxing the last set of 
constraints. A (M 3 +l) -dimensional multiplier vector (i.e., 



10 



*S A A 3 3 

<X> 3 (u 3 ) ^Minimize <t> 3 (z 3 ;u 3 ) = z2 12 z2 c i i i z i i i 

i 1 =0 i 2 = 0 i 3 =0 12 3 

A 3 [ #K A 3 1 



"12 3 12 3 



J 3 =0 3 



. i 1 =0 i 2 = 0 



L l J '2 J -3 



f A f 3 3l 3 A 2 

Minimize 2. 2. 2. ^iVa + u i 3 z UUi~H u i 

i,=0 i, = 0 i,=0 L 123 3-" 123 i =0 



15 



A 3 
Subject to * ' 3 

A £ z i,i,i, = l, i+l = l,..., 

i,=0 i,=0__ l\2,i 0 „^,f^ 



Uo s O 



, . 3 3 3 Ui = ° V° +1 V associated with 2 the 3 [th ]J£ 
u = \ u 0' u l' -' u Mj u 6 K 



20 constraint set will be denoted by with being fixed for 
notational convenience. (The zero multiplier is used for 
notational convenience.) Now, the three [ ] ^.dimensional 
assignment problem ([5] 6.1) is relaxed to a two-dimensional 
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assignment problem by incorporating last set of constraints 
into the objective function via the Lagrangian. Although any 
constraint set can be relaxed, we choose the last set of 
constraints for convenience. The relaxed problem is 



maximization of 0 3 (/"u 31 uf) with respect to the multiplier 
vector u 3 . It turns out that <P 3 is a concave, continuous, 

10 and piecewise affine function of the multiplier vector u 3 , so 
that the maximization of & 3 is a problem of nonsmooth 
optimization. Since many of these algorithms require a 
function value and a subgradient of & 3 at any required 
multiplier value (u 3 ) , we address this problem in the next 

15 subsection. We note, however, that there are other ways to 
maximize * 3 and the next subsection addresses but one such 
method . 




Subject to 




One of the major steps in the algorithm is the 
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VI. 2. Properties Lagrangian Relaxed Assignment 

w^wwww^v - 

Problem J" .1 

For a function evaluation of <P->. we show that an optimal 
(or suboptimal) solution of this relaxed problem ([Sj^.2) can 
be constructed from that of a two-dimensional assignment 
problem. Then, the nonsmooth characteristics of & 3 are 
addressed, followed by a method for computing a subgradient . 

i 2 ) an index 

• ' xri 3j 

10 c]; j_ J 3 = j 3 (i lf i 2 ) and a new cost function 
12 2 3 

c i 1 i 2 = c j 3 for i 2 ) * Joy: 

( [5] 6.3) 



Evaluation of A. 3 Define for each j^=~L(i^f 
j^j 3 (i lf i 2 ) ^argiHinjc^^^ + u^n i 3 = 0, 1, ...,My 



^V^iy^uj, fcr ( 3 i lf i 2 ) * (0,0), 
c o 2 o = £ Minimum \ 0 , c 0 3 0i + u? - 

i 3 =0 1 3 3; 

^Given an index pair (i lf i 2 ) , j 3 need not be unique, 
resulting in the potential generation of several subgradients 
15 ([5] 6. 11). (We further discuss this issue at the end of the 
Subsection [ 5.2.3].) 
Now, define 

[ .(5.4)] 
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i 1 i 2/ i 3 = 1 some (i^i,,!,) ;^ (u 3, =Minimize £ (2 2 ;u 3 } 2 £ £ 

' £^1, = !^ ^ = 1,...,^, (6.4) 

f> 2 

1^0 1 2 

z 2 ? ^6(0,1} for all i lf i 2 . 



As an aside, two observations are in order. The first is 
that the search procedure needed for the computation of the 
relaxed cost coefficients in ([5^.3) is the most 
computationally intensive part of the entire relaxation 
5 algorithm. The second is that a feasible solution [z 2] zi of a 
sparse problem ([5]£.l) yields a feasible solution /"z 2 ^zf of 
([5] 6^4) via the construction 

2 

rj_ ± if i 3 = J 3 and (i lf i 2 ) * (0, 0) ; 

) if i 3 *j 3 and i 2 ) * (0,0) ; 

L if c 3 Q0i ^u^0; 

.^3 3 rt 2 \l, if zf * ± =1 for some (i 

1 if c 00i + U; >0.z; i = < ^2^3 

00:L 3 2 3 ^ | 0 , otherwise. 
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Thus the two-dimensional assignment problem ([Sj^.4) has 
feasible solutions other than those with exactly one nonzero 
index if the original problem ([516^1) does. 

The following Theorem [53 6^1 states that an optimal 
5 solution of ([5] 6.2) can be computed from that of ([5] 6. 4). 
Furthermore, if the solution of either of these two problems 
is e-optimal, then so is the other. The converse is 
contained in Theorem [5] 6^2. 

Theorem [5] 6^1. Let [or 2 ] jbe a feasible solution to the 
10 two- dimensional assignment problem ([5] 6.4) and define w 3 
by 

32 ([5] 6.5) 

Wi l i 2 i, = w il i a if i 3 = j 3 and (i lf i 2 )* (0,0) , 

w li 2 i 3 = 0 if i 3 *j 3 and (i lf i 2 ) * (0, 0) , 

w ooi 3 =1 if c 0 3 oi 3 + Ui 3 3^0, 
"ooi 3 =0 if Co 0i3 + ul>0. 

Then w 3 is a feasible solution of the Lagrangian relaxed 
15 problem JJ^^gJ^and [03] 

[(w 3 ,^) = (p3 (w 2 ;u 3 )]<b 3 <p 3 {w 3 ' f u 3 ) =Q 3 = <fi 3 {w 2 ;u 3 ) . If, in addition, 
w 2 is optimal for the two [ ] ^dimensional problem, then w 3 is 
an optimal solution of the relaxed problem and 
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"1,1=2^ ^ii,!, for U^) * (0,0) 



12 i 3=0 1 2^3 



2 „ _ , . . v f 3 3 

L 3 ± 3 



^ 00 = 0 if (i x ,i 2 ) = (00) and Cqq^ + u.^0 for all i 3 



With ^the exception^ of one equality being converted to 

tfAA = l if (iir.i 2 ) = (00) a/idc Dni + uv *0 for some i .& {u 3 ) = S Au 3 ) 
an u inequality, the _followrng Jjttheorem is a_ converse of 



5 Theorem [5] 6.1. 

Theorem [5] 6.2. Let [w 2 ] w 3 be a feasible solution to problem 

^^^^ ^3— 

([5]£.2) and define w 2 by 

* ([516. [7] 6) 

- 3 (^ 3 ;u 3 )^<I>3(^ 2 ;u 3 )^ i2 =X; w z Wi for (i lf i 2 )M0,0) 



^ 0 2 0 = 0 if i 2 ) = (0,0) and c 0 3 0i3 + > 0 for all i 3 , 
r 0 2 0 = l if (i 1 ,i 2 ) = (0,0) and c 0 3 0i3 + u? 3 <: 0 for some i 3 . 



10 Then w 2 is a feasible solution of the problem ([5] 6. 4) and 
® 3 <p 3 (w 3 ;u 3 ) =<& 3 > <p 3 (w 2 ;u 3 ) . If, in addition, w 3 is optimal for 
([5] £.2), then w 2 is an optimal solution of ([5] £.4), 
$ 3 <p 3 (w 3 ;u 3 ) =® 3 = <p 3 {w 2 ;u 3 ) and & 3 (u 3 ) = S 3 {u 3 ) . 
[ The Nonsmooth Optimization Problem. ] 

15 

An Algorithm for Computing <P 2 and a Subgradient. 

[ ^^Given [u 31 uf the problem of computing ^ 3 ([u 3 ^ui) and a 
subgradient of ^([u 3 ^^) can be summarized as follows. 
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1. [ ^^^jStarting with problem ([5^.1), form the relaxed 

problem ( [5] 6 . 2) . 

2. To solve ([SJ^.2) optimally, define the ^two[ ] ^ 
dimensional assignment problem ([Sj^.4) via the 
transformation ( [5^6.3)^ 

3. Solve the two-dimensional problem ([5] 6.. 4) optimally. 

4. Reconstruct the optimal solution, say 
I 



w 

M, W- 



* 3 (u 3 )=<t> 3 (tf 3 ; U 3 ) = f; £ £ [Aii L i + £ \ u\ £ i: i -i 

i x = 0 i 2 =0 i 3 = 0 1 2 3 1 2 3 i 3 = o 3j. 1^0 i 2 = 0 12 3 

, of ([5] 6.2) via equation (6.5 [.6]) as in Theorem [5] 6.1. 
Then 



a 3 



M, M> 



0 3 (u 3 )^(tf 3 ;u 3 )=£ £ £ < V < V3 

ij=0 i 2 =0 i 3 =0 123 123 
A 3 f A A . 



15 6. A subgradient is given _by substituting [ iv 3 ] 



(6.7) 

[(5.10)] 
into 



the objective function ( [5] 6>2) , erasing the minimum sign, 

3 

and taking the partial with respect to u ± ^ . The result is 



1^ = 1^3 = 



d& 3 {u 3 ) 
duf 



£ £ 



i,=0 i, = 0 



for i 3 = 1, . 



(6.f 
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The next objective is to explain a recovery procedure, 
i.e., given a feasible (optimal or suboptimal) solution w 2 of 
([5] 6.4) (or w 3 of ([516.2) constructed in Theorem 

5 [5.L6.1), generate a feasible solution [z 3 1 z 3 of equation 
([516.1) which is close to w 2 in a ^sense to ^= ke specified. 
We first assume that no variables in ( [516. 1) are 

preassigned to zero; this assumption will be removed shortly. 
The difficulty with the solution w 3 in Theorem [5JJ5.1 is that 

10 it need not satisfy the third set of constraints in ([5jJL.l) . 
(Note, however, that if w 2 is an optimal solution for ( r 5 1 6. 4) 
and w 3 , as constructed in Theorem [5JJL-1, satisfies the relaxed 
constraints, then w 3 is optimal for ([5JJI.1).) The recovery 
procedure described here is designed to preserve the 0-1 

15 character of the solution w 2 of ( T5l 6 . 4) as far as possible: 
If and i^O or i 2 *0, the corresponding feasible solution z 3 of 
([516.1) is constructed so that for some (i 1 ,i 2 ,i 3 ). By this 
reasoning, variables of ^the^ form [,1 can be assigned a value 
of one in the recovery problem only if . However, 

20 variables will be treated differently in the recovery 
procedure in that they can be assigned ^either^ zero or one 

/ \ L 1~ 2 

independent of the value Ui 1 (l 2 ),i 2 (l 2 ) j, = qUW Q0 . This 



233 



NUMERI .07 USC 2 

increases the feasible set of the recovery problem, leading to 
a potentially better solution. 

Let be an enumeration of indices {i lf i 2 } of w 2 (or the 
first 2 indices of W ocristructed in equation (6^5) ) such that (w^ (I , ,i 2 (i 2 ))* ( o,oj = ^ 

5 for and w^ Q = lcalesymlOS (i x (1 2 ) , i 2 (I 2 ) )= (0, 0) for [I^OJI^O 

2 3 2 

regardless of whether Cj_ j_ ~Cj_ (i^) i 2 (l 2 ) i ^ or ^2 = ®* *' ^1^00 = 1 or 

not. (The latter can only improve the quality of the feasible 
solution. ) 

Next to define the two-dimensional assignment problem that 
10 restores feasibility, let 



2 2 

z2 c i ± z i ± 



i 2 =o i 3 =o 2^3 ^2^3 



A 2 

Subject to 2^f z l i = 1' 1 2 = 1, - , L 2 , 

i 3 =0 2 3 



£ Z2 i =1 ' V 1 '"*'^' C 2 2 i 3 =C i 1 (l 2 )i 2 (2 2 )i3 f ° r V 0 '~' L 1' 



i 2 =0 2 3 

e{0 f l} for all l 9 ,i.. 

([5] 6. [12] 9) 

15 Then the two-dimensional assignment problem that restores 
feasibility is 
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3 



£ t, a! 



( [5]g.l [4)0) 



iidj) i 2 U 2 ) i 3 



Minimize 





The next objective is to remove the assumption that all 

cost coefficients are defined and all zero-one variables are 

free to be assigned. We first note that the above construction 

of a reduced [dimension] t^^dime^ional assignment problem 

([5]6.1[3]1) is valid as long as all cost coefficients c 2 are 

defined and all zero-one variables in z 2 are free to be 

assigned. Modifications are necessary for sparse problems. If 

the cost coefficient is well-defined and the zero-one 

variable is free to be assigned to zero or one, then define 

as in ([5]6.[12]9) with being free to be assigned. 

Otherwise, is preassigned to zero or the corresponding arc is 

not allowed in the feasible set of arcs. To ensure that a 

feasible solution exists, we now only need ensure that the 

variables I 2 = 0, 1, I^z* 0 are free to be assigned for 
3 3 

z i 00' z 0i 0^2 ~ ® r 1' L i w i t]l the corresponding cost coefficients 
being well-defined. (Recall that all variables of the form 
and are free (to be assigned) and all coefficients of the 
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corresponding form and ^^M^^&^B^MS^^^' ^ I This is 
accomplished as follows: [ Il^f the cost coefficient is well- 
defined and is free, then define with being free. 

Otherwise, since all variables of the form and 

2 3 3 3 

5 c l Q = c i (l )00 + c 0i (I )O z oi 2 o are f ree to be assigned with the 

corresponding costs being well-defined, set ^ where the first 
term is omitted if i 1 (l 2 )=0 and the second, if i 2 (l 2 )=0. For 
i 1 (I 2 )=0 and i 2 (2 2 )=0 , define 



10 



15 



The recovery problem for the case [N=3]N = 3 is 



Minimize ). c l i z l i 

Y Y 2 3 ± 2 J -3 

A 2 

Subject to 2^ z l i = 1' i 2 = l,...,L 2 , 

i 3 = 0 2 3 

£ z li , = 1 ' 1 3 = 1 >~' M Z' 



2 2 =0 -"2^3 



2 2 3 

z 1 i 6 {0,1} for all 1 2 , i 3 .c 00 = c 000 



L 2-*-3 



20 



25 
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i=i (l),i 2 = i 2 {l 2 ) andY .=1; . . 
1 1 z * z z j. 2 i 3 Minimize 



erwise. 



2 

Ji 2 i 3 = l, 1 2 = 1, L 2 , 



i 2 =o 2 3 

z, 2 ± E {0,1} for all 2„i,. 

2 3 ^ J 
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<[5]6;1[4]1) 



Let Y be an optimal or feasible solution to this [2- 
dimensional] two-dimensional assignment problem. The recovered 
feasible solution z 3 is defined by 



_Jl, if i^i^I^ , i 2 = i 2 (I 2 ) and 
0, otherwise. 



( [5]6.1[5]2) 



^Said in another way, let Y|(I 2 (1 3 ) , i 3 (I 3 )} = 1^ =0 be an 
15 enumeration of indices {l 2/ i 3 } of Y such that 

3 . . J 1 ±f VV 2 12>' V^U*' VV^.y =1 fQr 

1 i- L 2 2 3 [o otherwise ' i 2 (i 3 > , i 3 d 3 ) 

(1 2 (1 3 ) t i 3 (l 2 ) I* (0,0) and (1 2 (1 3 ) , i 3 (1 3 ^) = (0 , 0) for 2 3 = 0 regardless 
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of whether Y 00 =l or not. Then the recovered feasible solution 
can be written as 

r 3 ^ 3 J 1 ' lf i l = i l (i 12)' i 2 = i 2( 1 l 2 )' i 3 = 1 3( i 3)' 

^ i i 2 2 i 3 \o, otherwise. 

( [5] 6.1 [6] 3) 



The upper bound on the feasible solution is given by 

10 

M, ML «, 

0 3 (u 3 ) ^<D 3 (u 3 ) <LV 3 V 3 (J 33 ) ZV 3 =J1 E E c/^*^ 



i 1= 0 i 2 =0 i 3 =0 



[V 3 ( 3 )=] (I5]6.1[7]4) 

15 

and the lower by <£ 3 ([u 3 ]ui) for any multiplier value ([u 3 ]uf) . 
In particular, we have 

I [(3)] 
where u 3 is any multiplier value, [(u 3 )]0^. is a maximizer of 

20 $ 3 ( [u 3 ] u*) , —3 is an optimal solution of the multidimensional 

assignment problem and 3 is the recovered feasible solution 

defined by ([5] 6. [16).] 

[ [Other Relaxations for the Multidimensional Assignment 
Problem. (ABP) ] 13^ 
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VII. OTHER RELAXATIONS FOR THE MULT I D I MENT I ONAL 
ASSIGNMENT PROBLEM 



In this section, we briefly describe other possible 
relaxations and their implications. These include linear 
5 programming relaxations and the corresponding duals. 

Recall from Section II that one starts with the definition 



of the zero-one variable z. . =1 and cost coefficients to form 



the problern 

[ a] 



10 



k . . . 
Zi Minimize J. c. . z. . 

* V^n 1 » 1 » 

Subject To z, . = 1 (i 1 = l / - f M L )r 

Va-N 

52 z i-i =1 (i 2 = l,-,M 2 ), . (7.1) 

■i V ... "J IN 



15 .... , X , vv 1 



(i p =l,-,M p and p = 2,-,N-l) 



f 



52 ^ ...i =1 (i w =lr- f M ) , 

j ii X l Z N W W 

- L l i 2" i N-l 

20 w h 

ere [c 0 ... 0 ]c 00 is arbitrarily defined to be zero. Here, each 
group of sums in the constraints represents the fact that each 
non-dummy report occurs exactly once in a "track of data". One 
can modify this formulation to include multiassignments of one, 
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some, or all the actual reports. The assignment problem 
([23^5) is changed accordingly. For example, if is to be 
assigned no more than, exactly, or no less than times, then 
the " = l" in the constraint ([2]^. 5) is changed to 

respectively . 
Modifications for group tracking 

_and multiresolution features of the surveillance region will 
be addressed in future work. In making these changes, one must 
pay careful attention to the independence assumptions that need 
10 not be valid in many applications. 

Again, the recent work of Poore and Drummond J^B^Poore 
§Si«i^JU««*Dr^^ 

multidimentional assignment problems. In D.W. Hearn, W.W. 



15 450^jy|gej^ 

E^V^^jsignif icantly extends the formulation of the track 
maintenance and initiation to new approaches. The discussions 
of this section apply equally to those formulations. 

2 0 In the linear programming relaxation, one replaces the 

zero-one variables constraints 

0<;z z <sl for all i x 6(0,1} for all i 2 , ...i^.., i w ( [6] 7.2) 



with the constraint 
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< I6]7. [2] 3) 

Then, the problem ([6]^7.1) can be formulated as a linear 
programming problem with the constraint ([6]^7.2) in ([6] 7.1) 
replaced by ([6] 7^.3) with a special block structure to which 
5 the Dantzing-Wolf e decomposition is applicable. Of course, 
after solving this problem, one must now recover the zero-one 
character of the variables in ([6]^7.1) and there are many ways 
to do this, such as using the two [ ] ^dimensional assignment 
problems. Commercial software is also available. 

10 

VI 1. 2. The Linear Programminq Dual and 
:ial Laorangian Relaxations, 



Given the linear programming relaxation, one can formulate 
the dual problem or the partial Lagrangian relaxation duals 

15 with respect to any number of constraints. In particular, this 
is precisely what is done in Section 3 on the Lagrangian 
relaxation algorithm presented. The much broader class of 
algorithms provided in the US Patent 5537119 (Aubrey B. Poore, 
Jr., Method and System for Tracking Multiple Regional Objects 

20 by Multidimensional Relaxation, US Patent Number 5537119, 
filed 14 March 1995, issue date of 16 July 1996) can be 
modified to remove the zero-one character when one relaxes M 
sets of constraints to an ^N-M[ dimensional] ^dimensional 
problem and recovers with an [N-M-j-1 dimensional] SM^Mf&J&X^? 
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dimensional problem. This avoids the problems associated with 
the NP-hard J^-AQ^ and [N-M+l] J^^+l]^ dimensional problems. 
However, to restore the zero-one character, one can do it 
sequentially with an assignment problem or with one of the many 
5 zero-one rounding techniques. These formulations are easy to 
work out and thus the details are omitted. 

The complete dual problem is another way of solving the LP 
problem and may indeed be more efficient in certain 
applications. In addition, the solution of this dual problem 

10 may provide excellent initial approximations to the multipliers 
for Lagrangian relaxations. 

[Multisensor Resolutions Problems and Their Solution via 
Linear Nonlinear Assignment Problems and Inequality 
Constraints. (ABP) The mathematical development for 

15 inequality constraints is understood by Aubrey B. Poore, but 
this work has not been written up. This feature is important 
in combining information from sensors with different 
resolutions such an IR and Radar.] 



2 0 [ Formulation and Algorithms for the Distributed Sensor 
Tracking: Track- to-Track Correlations. (ABP) 

Hot Starts for Track Maintenance and Initiation: Bundle 
Modifications (ABP) .] 
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YIJLLr, HOT STARTS FOR TRACK MAINTENANCE 

AND INITIATION: BUNDLE MODIFICATIONS 



Thus reuse of multipliers and the first proof that this 
reuse is actually valid was presented in this section. The 
5 reuse in track maintenance is demonstrated and discussed in the 
work of Poore and Drummond [xx] ^^^^^^^^^^^^^^^^^^^ 
Track^^^^ 

10 422^jtoston^^ The only 

item left is the modification of the bundle of subgradients for 
the use with the multiplier values as one goes through the 
recovery problem as in Section IV^ 3 . 6 or as one moves the 
window in track maintenance. It is this aspect of the 

15 nonsmooth optimization that adds an order of magnitude to the 
speed of the algorithms. This work is based on both a 
mathematical proof as in Section IV. 3. 6 as well as 
computational experience and heuristics. 

IX. [Adaptive Window Size Selection 

20 (ABP) -lADAPTIVE^^ 

These data association algorithms, based on the 
multidimensional assignment problem, range from sequential 
processing to many frames of data processed all at once! The 
data association problem for multisensor ^multitarget tracking 
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is formulated using a window sliding over the frames of data 
from the sensors. This discussion focuses on the work of Poore 
and Drummond and not on the formulations in Section [2] III. 
Firm data association decisions for existing track are made at, 
5 say frame k, with the most recent frame being k+M. Decisions 
after frame k are soft decisions. Reports not in the confirmed 
tracks are used to initiate tracks over frames numbered k-I to 
k+M. 

The length of these windows varies from sequential 
10 processing, which corresponds precisely to [^OJj^J^J^ an< ^ 

[M=1]M = 1, to user defined large values of I and M. [ . (In the 

^^^^^^^ 

case of sequential processing, we also have a temporary track 
file of dynamically feasible tracks, but incorrect data 
association.)] There is a marked change in performance over 

15 this range. As the two windows become exceedingly long, there 
is little statistically significant information gained and 
indeed performance can degenerate slightly. Thus, one must 
manually find the correct window lengths for performance in a 
given scenario and the algorithms do not change for a changing 

20 environment. Thus we propose an adaptive method for adjusting 
the window lengths. (The method has been highly successful for 
selecting the correct window length for multistep methods in 
ordinary differential equations.) 
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1. Compute the solution for different window lengths that 
overlap and differ by one or two frames. 

2. Compare the solution quality, i.e., the value of the 
5 objective function, for two adjacent window lengths. 

3. If there is a predefined improvement in either direction, 
we then, for stability, repeat the exercise for a shorter 
or longer on the first try. If there is consistent 

10 information, we adjust the window size in the indicated 

direction. This can be given a well defined mathematical 
formula in terms of the assignment problems of different 
dimensions . 

[• Algorithm Enhancements due to Data Structures (ABP) . 
15 • Construction and enumeration of the Best K-Solutions. 
• New Nonsmooth Optimization Techniques.] 

[ Sensitivity Analysis.] 

20 
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X. SENSITIVITY ANALYSIS 



For sequential processing in which the two [ ] ^dimensional 
assignment algorithm is used, one can use the LP sensitivity 
analysis theory and easily obtain the corresponding answers. 
5 For the multiframe processing, the optimal solution is not 
available; however, there are still two approaches one can use. 
First, the basic algorithm can perform the same sensitivity 
analysis at each stage (loop) of the algorithm as is done in 
the two [ ] ^dimensional assignment problem since the evaluation 
10 of function is equivalent to a two [ ] ^dimensional assignment 
problem. Alternately, one could use an LP relaxation of the 
assignment problem and base the sensitivity on the resulting LP 
problem. We currently see this as an important step in finding 
even better solutions to the assignment problem if so desired. 



15 [ 



Minimize 




Subject to 





20 



x i 6(0,1} 
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20^NLewAucy^^ EW AUCTION 



ALGORITHMS. 

In this section we present a new auction algorithm for the two 
^dimensional assignment problem. 
5 An important step in solving N-dimensional assignment 

problems for N >3, is finding the optimal solution of the [2- 
dimensional] tridimensional assignment problem. In particular, 
we wish to solve the [2 -dimensional] two^dime^ional problem 
which includes the zero-index. T [3 , ] gpically this problem can 

10 be thought of as trying to find either a minimum cost or 
maximum utility in assigning person to objects, tenants to 
apartments, or teachers to classrooms. We will follow the work 
of Bertsekas and call the first index set persons and the 
second index set objects. The statement of this problem is 

15 given below (11^1) when the number of persons is m and the 
number of objects is n. 

There are a couple of assumptions which we will make about 
20 (lJL^l) . First of all, we assume that lc i0 ]c± 0 and [c 0 j]c 0 j are 
well-defined and the corresponding variables [*io]gio anci 
[*oj]xoj are free to be assigned. Second if a cost c ± j happens 
to be undefined, then the corresponding variable will 
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be set to 0 . In effect if [c i j]c i j is undefined, we simply 
remove this cost and variable from inclusion in the problem. 

Notice that there are no constraints on the number of 
times person 0 and object 0 can be assigned. But notice that 
5 the first constraint requires that each non-zero person i must 
be assigned exactly once. Similarly, the second constraint 
forces each non-zero object j to be assigned exactly one time. 
Finally, the last constraint gives an integer solution, 
although we will see shortly that this constraint can be 
10 relaxed to admit any solution \Xa j>0] Xjj>0 . One reason for 
requiring that all of the costs of the form [c^lc^ and [c 0 j]c Oj - 
be defined is so that we are guaranteed a feasible solution 
exists for the given problem. 

RjsLa^cation^^ 

15 Relaxing the last constraint, we obtain the following 

problem: 

This lower bound is achieved and the second constraint in 
20 the original problem is satisfied if the following conditions 
are met : 
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For i=0, 



For i*0, 



1, if c 0J + Uj £ 0 
0, otherwise 



1, if j = arg min{c ik +u£\k = 0,l,...,h } 
0 otherwise 



All of j's are assigned only one time. 



The relaxation of the first constraint is analogous and would 
lead to similar results with i and j exchanged and the 
5 introduction of the multiplier u 1 . 

XI. 2. Relaxation of Both Constr^yjrtsT. ] 

This time we will relax both constraints and [using] use duality 

^^^^^ 

theory^to obtain a relationship between the multipliers u 1 and 
u 2 . 

10 [ a bunch of equations, theorems and the such here. 

From the above relaxation, we arrive at the following 
complementary slackness conditions which have been relaxed by 
a small value 6 . ] 

Definition: An assignment S and a multiplier set (u 2 ,u 2 ) are 
15 said to satisfy 6- complementary slackness (6 - CS) if 
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j L = arg min j c ik +u^\k eA{i) 



w ± = min 



{ c ik +u k\ k eA ( i > ' k * Ji) c ij + U i + Uj 2 = 0 for all (i,j)eS, 
. + Ui+u?^-e for all [i,j)6A. 



[ ] Forward Auction Algorithm. 

^^^^^ 

(1) Select any unassigned person i 

(2) Determine the following quantities: 

, iQ . caALIGNLj i = arg min |c i/c +u^|Jc c 
i I c ik +u%\k EA(i) , k* j.^. 

5 

In the selection of j ± above, if a tie occurs between 0 
and any non-zero index k f then select j ± as k. Otherwise, 
if there is a tie between two or more non-zero indices, 
the choice of j ± is arbitrary. 
10 Also if A(i) consists of only one element, then set w^ 00 . 

(3) Update the multipliers and the assignment: 
If j'i = 0, then 

(a) Add (i,0) to the assignment. 

(b) Update 
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If j ± * 0, then 

(a) Add (i/Ji) to the assignment. 

(b) Remove (i', Ji) from the assignment if j ± was 
previously assigned. 

(c) Update 




(d) Update 



10 Reverse Auction Algorithm. 

ij = arg min \ c jk + u k \k eB{j)f 

(1) Select^ any unassigned object j, such that 



15 



(2) Determine, the following* quantities : 
y.= min | c kj + u k \k €B(j), k* ij|c 0J + u) >0. 

:= -c^.i^arg min { c jJt + u* 1 1 /c 6B(j 



man 



In the selection of ij [ . ] above, if a tie occurs between 
0 and any non-zero index k, then select j ± as k. 
Otherwise, if there is a tie between two or more non-zero 
indices, the choice of j ± is arbitrary. 



251 




NUMERI .01USC2 

Also if B(j ) consists of only one element, then set 
Yj = 00 • 

(3) Update the multipliers and the assignment: 
If ij = 0, then 
5 (a) Add (0,j) to the assignment, 

(b) Update 
If ij * 0, then 



(a) Add to the assignment. 

10 (b) Remove (i^/jO from the assignment if was 

previously assigned. 

(c) U p d a t e 
V = '(%j^l)= -Y,-e^: = u^(Y,-P,)^= Y .-c^ + E f 

(d) U p d a t e 
15 If c oj + <0, then set u^:= -c oj u*:= -(c.^ + u^)= -y.-e . 

Combined Forward/Reverse Auction Algorithm 

1. Assume that u 2 is given as an arbitrary multiplier. 

2. Adjust the value of u 2 for each object j as follows: 

20 3. Run iterations of the Forward Auction Algorithm until 

all persons become assigned. 
4. Run iterations of the Reverse Auction Algorithm until 
all of the objects become assigned. 
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Note at the completion of the Forward auction step we have 
the following conditions satisfied: 

• for all objects j . 

• for all S. 

Thus we can prove the following proposition. 
Proposition: If we assume that at the start of the Forward 
Auction Algorithm and all of the persons are assigned via a 
forward step, then we have: 

for all A. 

c ij + u ) * ^ { c i* + 4} + 6 + u i l + u / = ° f ° ra11 

(i,j)e s. 

C.. + U* < min |c 1Jt +Uj?J + e for all (i,j)e S. 

Optimality of the Algorithm 

Theorem: 6-CS preserved during every forward and reverse 
iteration. 

Theorem: If a feasible solution exists, then the resulting 
solution is with mE of being optimal for the Combined 
Forward Reverse Algorithm. 



Implementation Specifics] 
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[Parallelization . ] 

[Here are but a few comments .] ^Z^^^^^^^^B^^^^SM, 
Although the algorithm appears to be serial in nature, its 
primary computational requirements are almost entirely 
5 parallelizable . Thus parallelization is planned. 

Step 2 is the computational part of the algorithm. 
Evaluating <& N _ k+1 and computing a subgradient use in the search 
procedure requires 99% of the computing time in the algorithm. 
This part uses two [ ] ^dimensional assignment algorithms, a 
10 search over a large number of indices, and a nonsmooth 
optimization algorithm. It is the second part (the search) 
that consumes 99% of the computational time and this is almost 
entirely parallelizable. Indeed, there are two [ ] ^dimensional 
assignment solvers that are highly parallelizable. Thus, we 
15 need but parallelize the nonsmooth optimization solver to have 
a reasonably complete parallelization. 

If a sensitivity analysis is desired or if one is 
interested in computing several near-optimal solutions, a 
parallel processor with a few powerful processors and good 
2 0 communication such as on the Intel Paragon would be most 
beneficial . 

The foregoing discussion of the invention has been 
presented for purposes of illustration and description. 
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Further, the description is not intended to limit the invention 
to the form disclosed herein. Consequently, variation and 
modification commensurate with the above teachings, within the 
skill and knowledge of the relevant art, are within the scope 
of the present invention. The embodiment described hereinabove 
is further intended to explain the best mode presently known of 
practicing the invention and to enable others skilled in the 
art to utilize the invention as such, or in other embodiments, 
and with the various modifications required by their particular 
application or use of the invention. It is intended that the 
appended claims be construed to include alternative embodiments 
of the invention to the extent permitted by the prior art. 
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TWhat is claimed is : ] 

[ A method for tracking a plurality of objects, 

comprising : ] 

[ repeatedly scanning a region containing a set consisting 
5 of one or more moving objects and generating N sequential 
images or data sets of said region, a plurality of observations 
in said images or data sets providing positional information 
for objects in said set;] 

[ determining a plurality of tracks, at least one track for 
10 each object in said set;] 

[ determining a plurality of costs, wherein each cost is for 
assigning one of said observations to one of said tracks ; 
defining a linear programming problem: 
Minimize ] 
15 [ Subject To (i 1 = 1,...,!^) 

(i 2 = 1, . . . ,M 2 ) 

(i p = 1, . . . ,Mp and p=2 , . . .N - 1) 
(i N = 1, . . . ,M N ) 

20 0 <; z u ... iN <l l for all ^,.,.4^] 

[wherein each c u _ iN is included in said plurality of costs, 
each i=l,...,N, being one of: (a) a number of observations 
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in an i th image or data set of said N sequential images or data 
sets; (b) a sum of a number of tracks in said plurality of 
tracks, and a number of said observations in the i th image or 
data set not assigned to one of said tracks; and (c) a number 
5 of tracks in said plurality of tracks;] 

[ solving said linear programming problem for values of 
2 u...iN for each i^.-ij,;] 

[ determining a value z il><#iN in {0,1} for each ii...^ 
corresponding to each z n iN , wherein said values z il#iiiN provide 
10 an optimal or near optimal solution to said linear programming 
problem; ] 

[ taking one or more of the following actions based on said 
optimal or near-optimal assignment of said plurality of points 
to said plurality of tracks:] 
15 [ sending a warning to aircraft or a ground or sea 

facility, ] 

[ controlling air traffic, 

controlling anti-aircraft or anti-missile equipment, 
taking evasive action, 
2 0 working on one of said one or more objects, 

surveilling one of said one or more objects.] 
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